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ABSTRACT 


We show how to formulate many continuous time-and-space search problems as gen¬ 
eralized optimal control problems, where multiple searchers look for multiple targets. 
Specifically, we formulate problems in which we minimize the probability that all of 
the searchers fail to detect any of the targets during the planning horizon, and prob¬ 
lems in which we maximize the expected number of targets detected. We construct 
discretization schemes to solve these continuous time-and-space problems, and prove 
that they are consistent approximations. Consistency ensures that global minimiz- 
ers, local minimizers, and stationary points of the discretized problems converge to 
global minimizers, local minimizers, and stationary points, respectively, of the orig¬ 
inal problems. We also investigate the rate of convergence of algorithms based on 
discretization schemes as a computing budget tends to infinity. We provide numer¬ 
ical results to show that our discretization schemes are computationally tractable, 
including examples with three searchers and ten targets. We develop three heuris¬ 
tics for real-time search planning, one based on our discretization schemes, and two 
based on polynomial fitting methods, and compare the three methods to determine 
which solution technique would be best suited for use onboard unmanned platforms 
for automatic route generation for search missions. 
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EXECUTIVE SUMMARY 


We develop formulations for continuous time-and-space search problems as general¬ 
ized optimal control problems, where searchers look for non-evading targets that move 
in a u'-dimensional area of interest, where w is a positive integer. We consider a large 
class of targets, where we assume the targets follow deterministic trajectories, given 
some information about their initial states or other parameters. We deal with two 
categories of targets, one where the targets coordinate their actions and one where 
the targets operate independently. 

We refer to the generalized optimal control problem when the targets coordi¬ 
nate their actions as the general target problem. We provide two different problem 
formulations for this category. The unconstrained and constrained problems ( GTP ) 
and ( GTP C ), respectively, are formulated to minimize the probability that all of the 
searchers fail to detect any of the targets during the planning horizon. The uncon¬ 
strained and constrained problems ( GTP e ) and ( GTP c,e ), respectively, are formulated 
to maximize the expected number of targets detected during the planning horizon. 

We refer to the generalized optimal control problem when the targets operate 
independently as the independent target problem, and formuate the unconstrained 
and constrained problems ( ITP P ) and ( ITP C ' P ), respectively, to minimize the prob¬ 
ability that all of the searchers fail to detect any of the targets during the planning 
horizon. We also formuate the unconstrained and constrained problems ( ITP e ) and 
(ITP c ' e ), respectively, to maximize the expected number of targets detected during 
the planning horizon. We develop discretization schemes to solve (GTP), ( GTP C ), 
(< GTP e ), (GTP c ' e ), (. ITP p ), (ITP C,P ), (. ITP e ), and (ITP c ’ e ), and show that the re¬ 
sulting finite-dimensional problems are consistent approximations to their infinite 
dimensional counterparts. Consistency of approximation ensures that global mini- 
rnizers, local minimizers, and stationary points of the discretized problems converge 
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to global minimizers, local minimizers, and stationary points, respectively, of the 
original problems. 

We consider a broad class of infinte dimensional optimization problems, which 
includes ( GTP ), ( GTP C ), ( GTP e ), (GTP c ’ e ), ( ITP *), (ITP c ’ p ), (ITP e ), and ( ITP c ' e ). 
We derive rate of convergence results for discretization methods used to solve prob¬ 
lems of this class by expressing the rate of convergence in terms of computational 
work. We find an upper bound for the rate of convergence by considering a finitely 
convergent optimization algorithm. The finitely convergent algorithm has no opti¬ 
mization error after a sufficiently large number of iterations, so the upper bound on 
the rate of convergence represents the best possible rate. We show that both super- 
linear and linear optimization algorithms are also able to attain this upper bound, 
and we identify specific discretization policies that achieve this best possible rate. 

We provide implementablc algorithms based on our discretization schemes 
along with numerical results to show that the problems ( GTP C ), ( ITP C,P ), and 
(JTP c,e ) are computationally tractable. For the 90 cases we consider, which include 
examples with three searchers and ten targets, we are able to compute near-stationary 
solutions in the range of 3 to 20 hours. 

We develop three heuristics for real-time search planning, one based on our 
discretization schemes, and two based on polynomial fitting methods, and compare 
the three methods to determine which solution technique would be best suited for use 
onboard unmanned platforms. Our numerical results indicate that our fixed-precision 
discretization scheme consistently provides solutions with the best objective values 
in the range of 40 to 200 seconds, and is therefore the best candidate for use as a 
real-time search planning method. 


xvi 



I. INTRODUCTION 

A. MOTIVATION AND BACKGROUND 

The need to protect a High Value Unit (HVU), such as an aircraft carrier 
or other large capital ship, from small boat attack is an important problem faced 
by navies around the world. The need to protect HVUs from these types of attack 
was brought into sharp focus following the attack on the USS Cole (DDG 67) in 
2000. The need was further demonstrated during a U.S. Navy war game in 2002, 
in which adversarial forces used swarms of small boats and aircraft to overwhelm a 
U.S. invasion fleet (Kahwaji, 2006). There are also indications that smaller countries 
could resort to an asymmetric strategy involving small boats if they were involved 
in a naval conflict against a major naval power (Kahwaji, 2006). A key component 
of any defense against small boat attack involves the search for and detection of 
potential adversaries. To quote Kahwaji, “Surveillance is key: If the raiders can be 
tracked as they swarm from their bases, they can be sunk with Rockeye cluster bombs 
and other munitions” (Kahwaji, 2006). In order to provide the best possible defense 
against small boat attack, it is clear that optimal utilization of search assets is highly 
desirable. 

Because the searchers’ trajectories have a significant impact on the probability 
of finding the targets within a given time horizon, we would like a way to determine 
the “best” trajectory for each of the searchers. This is a fundamentally difficult 
problem due to the nonlinearity and nonconvexity introduced by using probability of 
detection as the basic performance measure of a search platform. It becomes even 
more difficult if we consider multiple searchers looking for multiple targets. Given the 
increasing use of unmanned systems, such as unmanned aircraft systems (UASs) and 
unmanned surface vehicles (USVs), by militaries around the world, it is clear that 
there is a need for an automated method for finding optimal search trajectories given 
some initial intelligence information about potential adversaries. This need serves as 
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motivation for our work. In this dissertation, we consider several search trajectory 
optimization problems, including cases with multiple searchers looking for multiple 
targets. The searchers could be aircraft or surface vessels with sensors designed to 
detect targets that pose a threat to a HVU. We formulate these search problems as 
generalized optimal control problems, in continuous time and space. We develop and 
analyze algorithms based on discretization schemes to solve these problems; these 
schemes are consistent approximations in the sense of Polak (1997), and the resulting 
algorithms are guaranteed to converge to stationary solutions. 

B. SCOPE OF DISSERTATION 

We consider search trajectory optimization problems where searchers look for 
nonevading targets that move in a w-dimensional area of interest, where w is a positive 
integer. We limit the scope to searchers that move in continuous time and space, 
according to dynamics defined by ordinary differential equations. It is worth noting 
that this includes a large class of searchers, such as many manned and unmanned 
aircraft as well as many manned and unmanned surface vessels. We consider a large 
class of targets whose trajectories may not be continuous in time and space. While we 
only consider non-evading targets, targets are allowed to move intelligently. We allow 
for the case of targets that have perfect knowledge of the HVU’s position at all time, 
who then determine the trajectories required to strike the HVU in minimum time. We 
require that the targets’ motion is conditionally deterministic, which means that the 
targets follow deterministic trajectories given realizations of random variables that 
specify information about their initial states or other parameters. We assume that 
the probability distribution for each random variable is known, and we consider two 
situations. The first is when the random variables are dependent. This would be true, 
for example, if the targets chose to attack cooperatively in a swarm configuration. 
The second is when the random variables describing one target are independent of 
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the random variables describing another target. This would represent the case when 
the targets attack with no level of coordination, other than the desire to strike the 
same HVU. 

We refer to the generalized optimal control problem when the targets coordi¬ 
nate their actions as the general target problem. We provide two different problem 
formulations for the case when the targets coordinate their actions. The first formu¬ 
lation has corresponding problems ( GTP) and ( GTP C ), for the unconstrained and 
control-constrained cases, respectively, whose solutions minimize the probability that 
all of the searchers fail to detect any of the targets during the planning horizon. The 
second formulation has corresponding problems ( GTP e ) and ( GTP c,e ), for the un¬ 
constrained and control-constrained cases, respectively, whose solutions maximize the 
expected number of targets detected during the planning horizon. 

We refer to the generalized optimal control problem when the random variables 
are independent across targets as the independent target problem. We also provide 
two different problem formulations for the independent target category. The first 
formulation has corresponding problems ( ITP P ) and ( ITP C ' P ), for the unconstrained 
and control-constrained cases, respectively, whose solutions minimize the probability 
that all of the searchers fail to detect any of the targets during the planning horizon. 
The second formulation has corresponding problems ( ITP e ) and ( ITP c,e ), for the 
unconstrained and control-constrained cases, respectively, whose solutions maximize 
the expected number of targets detected during the planning horizon. 

Previous work on solving search problems using optimal control has focused 
on deriving necessary conditions for optimality (see, for example, Heilman, 1970, 
1971, 1972, and Lukka, 1977) in the tradition of Pontryagin, or sufficient conditions 
for optimality (see, for example, Hibey, 1982 and Ohsumi, 1991) in the tradition of 
Hamilton, Jacobi, and Bellman. Because it is unclear how to use these approaches to 
develop consistent approximations that converge to solutions of the original problem, 
we adopt the method developed by E. Polak (see, for example, Section 4.2 in Polak, 
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1997) to define optimality conditions. As in Section 4.3 of Polak (1997), we use Euler’s 
method to approximately solve the time-discretized ordinary differential equations 
governing the searcher dynamics. While it is possible to extend Polak’s method 
to include Runge-Kutta integration methods (see, for example, Schwartz & Polak, 
1996), doing so introduces additional complications in both theory and numerical 
implementation. Because this is the first attempt to utilize Polak’s method to derive 
optimality conditions and solutions for these types of generalized optimal control 
problems, we use Euler’s method throughout this dissertation to approximately solve 
the time-discretized ordinary differential equations. 

The problems we have formulated require spatial integration in addition to in¬ 
tegration in time. Spatial integration often requires more than one dimension, making 
numerical methods used to approximate spatial integrals more computationally ex¬ 
pensive than those related to time integration. We use Simpson’s integration rule 
as a higher-order approximation to the spatial integrals due to the fact that it helps 
limit approximation error, and is relatively simple to handle in both theory as well 
as implementation. 

C. LITERATURE SURVEY 

Modern search theory traces its roots to the formation of the United States 
Navy’s Operations Research Group (ORG) during World War II. During the war, 
Dr. B. O. Koopman and his associates in the ORG worked on improving the way the 
United States Navy conducted search during anti-submarine warfare operations (Iida 
et ah, 2002). Koopman (1980) is an unclassified and updated version of Koopman 
(1946), which summarizes the work done by the ORG during the war. Koopman 
(1980), and indeed much of the literature prior to the 1970s, focuses on stationary 
targets (Benkoski et ah, 1991). Research literature on moving targets can be grouped 
into two categories as given in Benkoski et al. (1991): 
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(i) articles which address special types of target motion, described below, that are 
amenable to analysis and 

(ii) articles which build on S. Brown’s work on conditioning with stationary tar¬ 
gets (see Brown, 1980) by developing general necessary and sufficient optimality 
conditions for moving-target problems. 

In this dissertation, we deal with target motion that is amenable to analysis, so we 
refer the interested reader to Benkoski et al. (1991) for a comprehensive survey of the 
literature related to category (ii). 

With respect to category (i), there are two special types of target motion that 
appear in the literature. The first type deals with targets whose motion is Marko¬ 
vian in nature, which means that the target moves in a random manner that can be 
modeled by a Markov process. Benkoski et al. (1991) again provides a comprehen¬ 
sive literature survey for Markovian target motion studies (until 1990). More recent 
research (see, for example, Washburn, 1998; Lau, Huang, & Dissanayake, 2008, and 
Sato & Royset, 2010) focuses on the development of specialized branch-and-bound 
algorithms for hireling an optimal path for the single searcher. In addition, Dell et al. 
(1996) present an exact procedure (utilizing a branch-and-bound algorithm) as well 
as six heuristics (local search, expected detection heuristics, genetic algorithms, and 
moving time-horizon heuristic) for solving the multiple searcher problem. 

The second type of special target motion in the literature is when the tar¬ 
get’s motion is assumed to be conditionally deterministic. The term conditionally 
deterministic means that the target’s trajectory depends on random variables, and 
if the random variables are given, then the target’s position is known for all time. 
This type of target motion has been investigated by many different researchers. In 
the case of discrete time, Royset and Sato (2010) presents a convex mixed-integer 
nonlinear program (MINLP) formulation for a route optimization problem involving 
multiple searchers who seek to detect one or more moving targets. Royset and Sato 
(2010) propose two solution approaches for the MINLP based on linearizations, one 
of which involves using a cutting-plane method. 
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For the case of continuous time, Stone and Richardson (1974) and Stone (1977) 
derive necessary and sufficient conditions to optimally allocate search effort to maxi¬ 
mize the probability of detection during a given planning horizon, lida (1989) extends 
the work of Stone and Richardson to include the problem of finding the optimal search 
plan which minimizes the expected risk (the expected search cost minus the expected 
reward), lida (1989) also derives the closed form of the optimal search plan when the 
target moves straight from a fixed point and selects its course and speed randomly. 
Pursiheimo (1976) derives a necessary condition for search plans to be optimal when 
the probability of detection is to be maximized and expected search time is to be 
minimized for a continuous time, discrete space model. The approach in Pursiheimo 
(1976) is noteworthy since it is the first to formulate the search problem for a target 
with conditionally deterministic motion as an optimal control model. The results in 
Pursiheimo (1976) are theoretical in nature, however, as no search plans are gener¬ 
ated. 

Following the work of Pursiheimo, there have been other formulations of the 
search problem using an optimal control model. Lukka (1977) uses an optimal con¬ 
trol model and derives a necessary condition for the optimal search plan assuming 
conditionally deterministic target motion. Mangel (1981) presents an approximate 
solution method, referred to as the “ray method,” that can be used in conjunction 
with the necessary conditions given by Lukka (1977) to determine the optimal search 
plan. Hibey (1982) and Ohsurni (1991) use an optimal control model to find sufficient 
conditions for the optimal search plan assuming Markovian target motion. Ohsumi 
(1991) also gives a method that can be used to numerically approximate the optimal 
search trajectories, with simulation studies assuming Markovian target motion. 

It should be noted that Lukka (1977), Mangel (1981), Hibey (1982), and 
Ohsumi (1991) all use an indirect approach to solving the optimal control prob¬ 
lem, employing the calculus of variations to obtain first-order optimality conditions, 
which result in a boundary-value problem that must be solved. A direct method 
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can also be used to solve the optimal control problem, which eliminates the need to 
solve a boundary-value problem. In a direct method, a time-discretization scheme is 
introduced, and the control problem is then transcribed into a nonlinear optimization 
problem that can be solved using standard techniques. Wasburn (1990) employs a 
direct method to construct a tactical decision program, JITTER, which is designed to 
find an optimized trajectory for a submarine transiting from a given starting location 
to a given terminal location, while minimizing the total acoustical energy received by 
n listeners. JITTER utilizes discretization and a steepest-descent method that grad¬ 
ually modifies a user-supplied initial trajectory into an optimized path by making 
first-order corrections. 

A more recent approach, given in Chung et al. (2010), uses a direct method 
based on Chapter 4 of Polak (1997) to solve an optimal control problem which Ends 
optimized search plans for a target that moves straight down a channel at constant 
speed. The results in Chung et al. (2010) indicate that the direct method based on 
Chapter 4 of Polak (1997) can be used to generate optimized search plans for as many 
as three searchers. Because the direct method based on Chapter 4 of Polak (1997) 
demonstrates the ability to produce search plans and allows us to develop consistent 
approximations, we focus on this method in the dissertation. 

D. CONTRIBUTIONS 

In this dissertation, we extend the work of Chung et al. (2010) in many ways. 
Our target motion model is more general than that found in Chung et al. (2010), and 
as a result the framework we develop allows for the formulation of many important 
continuous time-and-space search problems as generalized optimal control problems. 
We provide four different examples of problem formulations, namely ( GTP ), ( GTP e ), 
( ITP P ), and (JTP e ), that can be modeled using our framework. We also develop a 
min-max formulation that can be used in future studies to find searcher trajectories 
that minimize the maximum expected time of first detection of all the potential 
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targets, as well as a model that can extend our detection-based approach to one that 
includes herding potential threats away from a high value unit. As in Chung et al. 
(2010), we also use a direct method to solve the generalized optimal control problems 
based on Chapter 4 of Polak (1997), but we are the first to provide proofs to show that 
the discretization schemes are consistent approximations to their infinite-dimensional 
counterparts. 

We develop rate of convergence results based on expressing the rate of conver¬ 
gence in terms of computational work rather than the traditional number of iterations 
or level of discretization. Using this approach, we provide an upper bound on the 
rate of convergence that can be achieved by any optimization algorithm. We also 
provide discretization policies for superlinearly and linearly convergent optimization 
algorithms that achieve this upper bound. In addition, we use the rate of convergence 
results we obtain to provide insight regarding the choice of numerical method used 
to approximately solve the differential equations as well as approximate the spatial 
integration when solving generalized optimal control problems. 

We provide implementable algorithms based on our discretization schemes, 
and use them to produce numerical results to show that our method is computation¬ 
ally tractable. While Chung et al. (2010) considers a single target and as many as 
three searchers, we give numerical examples that include three searchers looking for 
as many as ten targets. We also investigate methods to reduce the computational 
cost necessary to obtain these numerical solutions, including adaptive discretization 
schemes and heuristics based on polynomial fitting. 

E. ORGANIZATION 

The remainder of this dissertation is outlined as follows. Chapter II gives the 
objective function definitions, describes the target motion model, defines the searcher 
dynamics, and provides preliminary generalized optimal control problem formulations. 
Chapter III introduces the spaces necessary to complete the problem formulations and 



develops consistent approximations for the generalized optimal control problems we 
consider in the dissertation. In Chapter IV, we develop rate-of-convergence results 
for different classes of optimization algorithms that can be used to solve generalized 
optimal control problems. Chapter V provides implementable algorithms and gives 
our numerical results. In Chapter VI we present conclusions and suggest future 
research opportunities. The Appendix provides some mathematical background. 
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II. SITUATIONAL DESCRIPTION AND 
PROBLEM FORMULATIONS 

We consider the situation depicted in Figure 1, where a High Value Unit 
(HVU) is operating in a two-dimensional area of interest. The HVU follows a fixed 
trajectory, {x°(t) G M 2 : 0 <t <T < oo}, during the finite planning horizon [0, T\. 
Without loss of generality, we assume throughout the dissertation that T = 1, and 
therefore our normalized planning horizon is [0,1]. During the transit from a;°(0) 
to x°(l), the HVU is under threat of attack from L targets. We assume that the 
motion of these targets is conditionally deterministic, i.e., the I th target’s position, 
y l (t;a ), is known for all time t e [0,1] given the realization of a random vector a 
which takes values in a compact set A C HU", where w is a positive integer, ft should 
be noted that a will be used to represent both a random vector and its realization, 
but the nature of a should be clear from the context. We assume that the probability 
distribution for a is known. In an effort to detect the potential threats to the HVU, 
there are K searchers operating in the vicinity of the HVU. Our goal is to determine 
the trajectory, x k (t), for each of the searchers to follow during the planning horizon 
that optimizes objective functions we define in Sections II.A and II.B. 

There are many different ways to measure the effectiveness of a proposed 
search plan. In this chapter, and throughout the dissertation, we focus on two types 
of objective functions. The theory we develop can be extended to include other types 
of objective functions (see, for example, Chapter VI, Sections B.l and B.2), but for 
our analysis and numerical results we concentrate on two types of objective functions: 
the first type seeks to minimize the probability that all of the searchers fail to detect 
any of the targets during the planning horizon; the second type seeks to maximize 
the expected number of distinct targets detected during the planning horizon. 

In this chapter, we begin by developing both types of objective functions for 
the case of statistically dependent targets. The term dependent targets refers to the 
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case when the random variables that define the targets’ motion are dependent among 
targets. A coordinated swarm attack against the HVU would be appropriately mod¬ 
eled by the dependent target case. We then derive both types of objective functions 
for the case of independent targets. The term independent targets refers to the case 
when the random variables that define the targets’ motion are independent across 
targets. It should be noted that in the case of independent targets there may be 
dependence between the random variables that specify information about the pa¬ 
rameters for a particular target. An attack by multiple targets whose only level of 
coordination is the desire to strike the same HVU would be appropriately modeled 
by the independent target case. Next, we discuss the target motion model. Finally, 
we give preliminary definitions of the problems we consider in this dissertation. We 
complete these problem definitions in Chapter III, after we define the appropriate 
spaces for the searcher control inputs. 
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A. COORDINATED TARGETS 

Given a G A C M 1 " common to all targets, the I th target follows a deterministic 
trajectory, {y l (t; a) G M ny : 0 < t < l}, where n y is a positive integer. As an example, 
a could be used to represent the uncertainty in target starting location and starting 
time. Then, given the starting location and starting time for the I th target, some 
type of deterministic algorithm can be used to find the trajectory for the I th target. 
While the formulation that follows is valid for any w, the dimension of A, it appears 
difficult to carry out the analysis for an arbitrary w. This is because w influences 
the choice and characteristics of the numerical integration scheme used to calculate 
search plans. During later proofs of convergence it is necessary to quantify the error 
introduced by the numerical integration scheme, so we now fix w = 2. 

We let r k ’ 1 : M™ x M ny —» [0,oo), were n is a positive integer, denote the 
detection rate for searcher k against target /, which is defined such that r k,l (x k ,y l )At 
approximates the probability of the k th searcher in state x k G M n detecting the I th 
target in state y l G M ny during a small time interval [t,t + At). The states for the 
searcher and target typically involve their locations, but could also include other 
quantities such as heading and time of day. The detection rate reflects the sensor 
effectiveness and we typically have that r k,l (x k : y l ) is some decreasing function in 
the “distance” between x k and y l . We put distance in quotations because when 
n ^ n y . it is necessary to omit certain portions of the state vectors for x k and y l in 
order to properly define a norm that provides a measure of the closeness of the two 
state vectors. The detection rate, r k,t (•,•), can be selected to reflect various types 
of sensors and their performance against different types of targets. For theoretical 
and computational reasons, r k,l ( •, •) must satisfy certain differentiability assumptions, 
which we state in Assumption 111.3. 

Next, we derive the detection model for a particular a, but first we define some 
notation. In a manner similar to Chung et al. (2010), given a particular trajectory 
for the k th searcher {x k {t) : 0 < t < 1} and a particular trajectory for the I th target 
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{y l (t\a) : 0 < t < 1} under realization a, we denote the probability that the k th 
searcher does not detect the I th target during [0, t], t 6 [0,1], by q k,l (t\ a). We assume 
that events of detection in non-overlapping time intervals are all independent, so we 
can calculate q k ’\t] a) recursively using the difference equation 1 

q k,l (t + At] a) = q k ' l (t;a ) (l - (r k,l (x k (t), y l (t) a)) At + o(At))) ,q k,l (0;a) = 1, (II. 1) 


which becomes the parameterized differential equation 


q k -‘(t ]a ) 4 a) = -yy ;a )r‘V(i),i/(i;a)), q ki ( 0;a) = 1, (02) 


as At —> 0, with solution 


q k '\t-a) = exp ( - / r k ’\x k (s), y\s- a))ds ) . 

Jo 


(H3) 


We assume that the searchers make independent detection attempts and can 
simultaneously detect multiple targets, and hence it follows from (II.3) that the con¬ 
ditional probability that no searcher detects any target during the time period [0,1], 
given a and collection of searcher trajectories, {x k {t) : 0 < t < 1}, k = 1, 2, is 

simply the product 


KL, 1 \ / K L .1 

11II ex P ( _ / r k ’ l (x k (t), y\t\ a))dt ) = ex P ( “ \ 

k =1 i=i \ Jo / \ k=1 l=1 Jo 

= ex p - / 

\ J° k =l l=i 


r k,l (x k (t),y l (t] a))dt 
r k, \x k {t ), y l (f; a))dt 


Similarly, the conditional probability that no searcher detects the I th target 
during the time period [0,1], given a and collection of searcher trajectories, { x k (t ) : 
0 < t < 1}, k = 1, 2, K, is given by the product 


K 

H exp 


k= 1 


"1 \ / r 1 K 

j\ „.l/ 


r k ’ L (x k (t),y L (t-,a))dtj = exp I — J ^^ r k '\x k (t),y l (t; a))dt J . 


1 Recall that if a function / : 
on page 304 in Ross (2007). 


k =1 


is o(x), then lim^^o = 0; see for example Definition 5.2 
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The conditional probability that at least one of the searchers detects the I th target 
during the planning horizon [0,1], given a and collection of searcher trajectories, 
{x k (t) : 0 < t < 1}, k — 1, 2,..., K, is given by 


r l K 

1 - exp ( - / r k ’\x k (t),y l (t- a))dt ) . 

Jo k =i 


(11.4) 


Then, the expected number of targets detected during the time period [0,1], given a 
and collection of searcher trajectories, {x k (t) : 0 < t < 1}, k — 1, 2,..., K, is given by 


E 

i=i 


»1 K 


1 - exp - / y ^r k ' l (x k (t),y l (t] a))dtt 


k =1 


( 11 , 5 ) 


We let (j) : A —> M be the probability density function of a. For theoretical and 
computational reasons, 0(-) must satisfy certain differentiability assumptions which 
we state in Assumption 111.1. Then, the probability that all of the searchers fail to 
detect any of the targets during [0,1] is given by 




1 I< L 

exp ( — / ) <p(a)da , 

Jo k =i i=i 


(IL6) 


and the expected number of targets detected during the time period [0,1] is given by 


E 


'" eA «=i 


„1 K 


1 — exp f — Y rk,l ( xk (^y l ^ a )) dt 


E 

(=1 


k =1 
1 K 


(j)(a)da 


1 - 


' a£A 


exp J r k, \x k (t ), y l (t ; a))dtj 4>(a)da 


( 11 . 7 ) 


B. INDEPENDENT TARGETS 

We now assume that the random variables that the target motion is con¬ 
ditioned upon are independent among targets. We then have a vector of random 
variables a 1 e A C M> w , one for every target, and we assume that the probabil¬ 
ity distribution for a 1 is known for all l = 1,2,...,L. Given a particular realiza¬ 
tion of the random vector, a 1 , then the I th target follows a deterministic trajectory, 
[y l (t] a 1 ) G W ly : 0 < t < l}. We note that the components of a 1 can be dependent. 
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Following the same development as (II.1) through (II.3), with a replaced by 
a 1 , we find that the conditional probability that the k th searcher fails to detect the 
I th target during the time period [0,1], given a 1 and searcher trajectory { x k {t ) : 0 < 
t < 1}, is given by 

exp r k ’ l (x k (t),y l (t-a l ))d?j . (II.8) 

We assume that the searchers make independent detection attempts, and hence it 
follows from (II.8) that the conditional probability that no searcher detects the I th 
target, given a 1 and collection of searcher trajectories, {x k {t) : 0 < t < 1}, k — 
1, 2,..., K , is simply the product 

TT exp f — f r k ’ l (x k (t),y\t-a l ))dt] = exp f - V f r k \x k {t),y\t-a l ))dt 

k =i V Jo J \ k=1 J o 

( r 1 K \ 

■ /o k= i 1 


= exp 


We let (j) 1 : A —> M be the probability density function of a 1 , l = 1, 2,..., L, and 
require it to satisfy certain differentiability assumptions, which we state in Assump¬ 
tion III.31. Then, the probability that all of the searchers fail to detect the I th target 
during the time period [0,1] is given by 

f y^r k ’ l (x k (t),y l (t-,a l ))dt \ <f) l (a l )da l . (II. 9 ) 

/o k =i J 

Hence, the probability that all of the searchers fail to detect any of the targets during 

[0,1] is given by 


/ a l eA 


exp 


L 

n 

i=i 


/ a‘eA 


exp 


— f y ^ y r k ’ l (x k (t),y l (t] a l ))dt \ <j) l (a l )da l 

k =i 1 


( 11 . 10 ) 


Based on (II.9), the probability that at least one of the searchers detects the I th target 
during the planning horizon [0,1] is given by 

' rl K \ 

— j yy r k ’ 1 (x k (t) , y l (t; a l ))dt J (f)\a l )da l . (II.11) 


exp 


I a l eA 
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Then, the expected number of targets detected during the time period [0,1] is given 
by 

exp - / 
eA y Jo 

We use (II.6), (II.7), (II.10), and (11.12) to formulate the generalized optimal 
control problems that we consider in this dissertation, but first we discuss the target 
motion model which we use to generate the target trajectories {y l {t\ •) G : 0 < 
t < 1}. 

C. TARGET MOTION MODEL 

We develop this section for the coordinated target case, but note that the inde¬ 
pendent target case is identical except that a is replaced by a 1 . Given a particular a, 
the target trajectories, {y l (t; a) G M ny : 0 < t < 1}, l = 1, 2, could be generated 
in numerous ways. In Chapter III, we make very light assumptions regarding ?/(•; •), 
so the theory we develop applies to a broad class of targets. In order to provide a 
specific numerical example, we make the conservative assumption that the targets 
have full knowledge of the position of the HVU at all time. Based on this information 
about the HVU, we model the targets as Dubins vehicles, i.e., nonholonomic vehicles 
that are constrained to move along planar paths of bounded curvature, without re¬ 
versing direction, that act intelligently and follow trajectories that seek to minimize 
the time required, tf, for each of them to hit the HVU. The motion of the targets is 
subject to additional constraints, which we now discuss in detail. 

For any t G [0, tf], let n y = 3 and y l (t;a ) = (y[(t; a), y ] 2 (t; a), y l 3 (t; a)) T G M 3 , 
where T denotes the transpose of a vector, be the state of the I th target at time t, with 
y[(t; «)gR and y l 2 (t', a)Gl denoting the horizontal and vertical components of the 
location of the I th target, respectively, and y l 3 (t',a ) G M denoting the heading of the 
I th target measured from the horizontal axis at time t. For any t G [0, tf], the control 
input, u l,tar (t), for the target is the rate of change of the heading, which is restricted 


K 

V 

k =1 


r k,l (x k (t),y l (t;a l ))dt ] 4>\a l )da l 


(II-12) 
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to the range [— u l,tar , u l,tar ]. The target’s speed, \\(y[(t-,a),y l 2 (t-,a))\\, is restricted to 
the range [v l min , v l max \. We require that the target’s initial and final speeds, v l 0 and v l f, 
be specified. Recall that a is based on a distribution, then given a particular a, the I th 
target’s starting position, y l 0 (a), and starting time, t l 0 (a), are known and we generate 
its trajectory by approximately solving the following optimal control problem: 

min t f 

s.t. y[(t-,a) = \\(y[(t;a),y l 2 (t-,a))\\ cosy l 3 (t; a), Vt G [t l 0 (a),t f ] (11.13) 

f&fa 01 ) = O'), y l 2 (t-,a))\\ sin y l 3 (t-, a), Vt e [t l 0 (a),t f ] (11.14) 


y l 3 (t;a) = u l ’ tar (t),Vte[t l 0 (a),t f } (11.15) 

V l min < <V l max , Vt G [4(«), tf] (11.16) 

-U l ’ tar < U l ’ tar (t) < U l ’ tar , Vt e[t l 0 (a),t f ] (11.17) 

v l o = ||(^(to(a);a),^(4(«);«))ll ( IL18 ) 

v l f = \\(y[(tf,a),y l 2 (tf,a))\\ (11.19) 

y l (t l 0 (a);a) = y l 0 (a ) ( 11 . 20 ) 

(xi(tf),x° 2 {t f )) = (y[(t f ;a),y l 2 {tf,a)) ( 11 . 21 ) 


The target dynamics are given by equations (11.13), (11.14), and (11.15). Constraint 
(11.16) restricts the target’s speed to the range [v l min , v l max \. Constraint (11.17) restricts 
the target’s turn rate to the range [—u l,tar ,u l,tar ]. The boundary conditions given by 
(11.18) and (11.19) require the target to start and end at the given initial and final 
speeds, v l 0 and v l j. The boundary condition (11.20) requires the target to start at 
the starting position, y l 0 (a), at the target’s starting time, t l 0 (a). The final boundary 
condition (11.21) requires the target to hit the HVU at the final time, tf. 

We solve this optimal control problem via a direct method that fits seventh- 
order polynomials to the target trajectories in a manner similar to that found in 
Yakimenko (2000) and Ghabchcloo et al. (2009). We note that the optimal control 
problem used to obtain the target trajectories is time-optimal, and we are concerned 
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with search during the planning horizon [0,1]. This will require truncation of the 
target trajectories after one hour of travel time, and is further discussed in Section 
V.B.l. 

We are now ready to state the generalized optimal control problems that we 
consider in this dissertation. 

D. GENERALIZED OPTIMAL CONTROL PROBLEMS 

Our goal is to find trajectories for the K searchers that optimize the expressions 
(II.6), (II.7), (II.10), and (11.12). We assume that the motion of the k th searcher, 
which we refer to as its “physical” dynamics, is governed by the differential equation 

x k (t) = h k (x k (t), u k (t)), t G [0,1], x fc (0) = (11.22) 

where £ k is the vector of initial conditions for the k th searcher, which could include 
things such as the initial position and the initial heading of the searcher, the control 
u k {t) G is the control input to the k th searcher at time t, which could be the rate 
of change of the heading of the searcher, and therefore h k : M n x M mfc —> M n . We also 
define m = Y2k =i m £ = (£ 1 > —■> £ A ) T , and u(t) = (u l (t) T ,..., u K ( t) T ) T , and therefore 
the collection of h k is given by h : M. nK x M m —> W nK . We can also handle time- 
varying systems and arbitrary planning horizons by using standard transcriptions; 
see for example p. 493 in Polak (1997). 

In order to completely state the generalized optimal control problems, it is 
necessary to limit the searcher controls to certain spaces which we define in Chap¬ 
ter III. For now, we give preliminary definitions of the generalized optimal control 
problems. 

We define x v,k (t ), t G [0,1], as the solution to (11.22) with r/ = (£,w) as the 
input, and let x v (t) = [x v,1 (t) T , ..., x n,K {t) T ) T ■ We assume that the solution x v,k (t), 
t G [0,1], exists and is unique based on assumptions we will make in Chapter III. For 
coordinated targets, we refer to the problem that minimizes the probability that all 
of the searchers fail to detect any of the targets during [0,1] by choosing the best rj 
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as the general target problem ( GTP ). From (II.6), we have the following generalized 
optimal control problem formulation: 


(GTP) : min 

v 



exp 


/ ,1 k L 

-/ ££<-‘v* 

\ • j/ ° fc=! ;=i 


x v,k (t),y (t] a))dt (j)(a)da > . (11.23) 


We note that constraints on searcher controls and initial conditions will be given in 
Chapter III. The constrained general target problem ( GTP C ) is exactly the same as 
(GTP), with the addition of 


u(t) eU,te [0, l], 


(11.24) 


where U is a convex, compact subset of M m . For coordinated targets, we refer to 

the problem that maximizes the expected number of targets detected during the time 

period [0,1] by choosing the best rj as the general target problem (GTP e ). From 

(II.7), we have the follwing generalized optimal control problem formulation: 

k 


(GTP e ) : min £ 


i=i 


exp 




f y ftr k,l (x v ’ k (t),y l (t-,a))dt) 4>(a)da 
k= i ' 


(11.25) 


The constrained general target problem (GT P c,e ) is exactly the same as (GTP) with 
the addition of (11.24). 

For independent targets, we refer to the problem that minimizes the probability 
that all of the searchers fail to detect any of the targets during [0,1] by choosing the 
best 7] as (ITP P ), and the problem that maximizes the expected number of targets 
detected during the time period [0,1] by choosing the best y as the independent 
target problem (ITP e ). From (II. 10) and (11.12), respectively, we get the following 
generalized optimal control problem formulations: 


(ITP P ) : min n 


. i=i 


,1 K 


x l £A 


exp 


r k ’ 1 (x v,k (t), y l (t; a l ))dt (ft(a 1 )da 1 


and 

(ITP e ) : min <j ^ 


i=i 


k=l 


- J ex p (- / X 

Ja l &A \ JO fc=1 


(11.26) 


r k ’ l (x v ' k (t),y l (t] a l ))dt (ft(a l )da l 


(I 


,27) 
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The constrained independent target problems ( ITP C,P ) and (PTP c,e ) are exactly the 
same as ( ITP P ) and (PTP e ), respectively, with the addition of (11.24). We again note 
that constraints on searcher controls and initial conditions for (PTP P ) and (JTP e ) 
will be given in Chapter III. 
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III. 


CONSISTENT APPROXIMATIONS 


The problems ( GTP ), ( GTP C ), ( GTP e ), ( GTP c ' e ), ( ITP P ), ( ITP C ' P ), (JTP e ), 
and ( ITP c,e ) defined in Chapter II are infinite dimensional in both time and space. 
In order to numerically solve these problems, some form of discretization is necessary. 
For the discretization scheme to be useful, it must lead to implcmentable algorithms 
that can solve the discretized problems in a reasonable amount of time, and the result¬ 
ing solutions must correspond to the solutions of the original problems in some sense. 
The discretization schemes provided in this chapter meet these requirements. In this 
chapter, we also formally define the term consistent approximation, used to describe 
the relationship between solutions obtained using finite-dimensional approximating 
problems to the solutions of the original infinite-dimensional problems. 

As alluded to in Chapter II, the problems (GTP), ( GTP C ), ( GTP e ), ( GTP c ’ e ), 
(. ITP P ), (ITP C)P ), ( ITP e ), and ( ITP c,e ) are well defined only if the allowable searcher 
controls are restricted to certain spaces. This chapter begins with specific definitions 
of those spaces. The remainder of the chapter is then split into two main sections, 
beginning with the treatment of coordinated targets and followed by independent tar¬ 
gets. Both of these main sections proceed in a similar manner, as follows. First, we 
define an “information state,” which we use to write the generalized optimal control 
problems in terms of the spaces defined in Section III.A. Next, we state our assump¬ 
tions and define optimality conditions for the generalized optimal control problems. 
Then, we develop consistent approximations for the time-discretized search problems. 
Finally, we show that the time- and space-discretized search problems are consistent 
approximations for the original, continuous time-and-space search problems. 

A. CONTROL INPUT 

We define the problems (GTP), (GTP C ), (ITP P ), (ITP C ’ P ), (ITP e ), and 
(ITP c,e ) using the following spaces, which are adopted from Chapters 4 and 5 of 
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Polak (1997) and described here for the sake of completeness. For the sake of brevity, 
we do not explicitly deal with the problems ( GTP e ) and ( GTP c,e ) in this chapter. 
We note, however, that results developed in Section III.C for the problems ( ITP e ) 
and ( ITP c,e ) can be trivially extended to include ( GTP e ) and ( GTP c,e ), respectively, 
if a 1 is replaced by a and 4>\-) is replaced by 0(-). We assume that the control u is 
an element of a subset of L™[ 0,1], the space of Lebesgue square-integrable functions 
from [0,1] into M m . Hence, initial condition and control input pairs are elements of 
the space 

H 2 = R n x L%[0,1]. (III.l) 

We denote the L2 n [0,l] inner product for any pair of functions u,v G L™[ 0,1] by 
(u,v ) 2 = fy (u(t),v(t)) dt, with (•, •) denoting the Euclidean inner product. The 
corresponding norms ||u|| 2 = ( u,u ) l J 2 and || • || — (-,-) 1 ^ 2 . For any 77 = (£,w) G H 2 , 
with £ G M n and u G L™[0,1], and any rj' = ( ') G H 2l with G and 
u' G L™[ 0,1], we define the inner product and norm on H- 2 , respectively, by 


(v,v')h 2 = it,?) + {u,u') 2 

(I 11 . 2 ) 

IWk = (M 2 + Wl2) 1/2 - 

(III.3) 


We further restrict the control u to be an element of L™[0,1], the space of 
essentially bounded, measurable functions from [0,1] into M m , but find it convenient 
to retain the inner product and norm of T™[0,1]. Hence, as in Section 5.6 of Polak 
(1997), we define the pre-Hilbert space 

H^ 2 ^R n xLZ, 2 [0,l], (III.4) 

where L™ 2 [0,1] denotes the space 0,1] equipped with (•, -) 2 and || • || 2 . We note 
that II.y- 2 is dense in H 2 . 
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We also assume that there exists a p max < oo such that controls of interest at 
all times are located in the interior of f?(0, p max ) — {v G M m | ||n|| < p max }. Therefore, 
we focus on the space 

H = IT x U C # 00,2 (III.5) 

of initial conditions and control inputs, where 

U 4 { u G Z£ )2 [0,1] | u(t) G #(0,p max ), Vt G [0,1]} . (III. 6 ) 

When dealing with differentiability statements we restrict ourselves to the 

subset 

H° = M n x U° C H. (III. 7) 

where 

U° - {«G L™ 2 [0,1] | u(t) G #(0, 7 /w), Vt G [0,1]} , (III. 8 ) 

and 7 G (0,1) is close to one. It is clear that, when p max —* oo, H° “fills ’ 7 H^- 

Finally, we allow for the inclusion of a specific type of pointwise control con¬ 
straint of the form u(t) G U for all t G [0,1], where U C M m is a convex compact 
subset of £>(0, 7 p max ). This means we only consider pointwise control constraints of 
the form u G U c , where 

u c = {« G U° I u(t) euc B(0, 7 Pmax), Vt G [0,1]} . (III.9) 

We also use the notation 

H c = PxU c . (III. 10) 

B. COORDINATED TARGETS 

1. Information State and Optimal Control Problems 

In this section, we define an “information state” which we use to write the 
generalized optimal control problems in terms of the spaces defined in Section III.A. 
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For a given rj £ H. recall that x v,k (t), t £ [0,1], is the solution to (11.22) with rj as the 
input. Based on (11.23) we define the objective function / : H —y M for any 77 £ H by 


K L 


f(v) = exp - EE r k ’\x v,k (t),y l (t;a))dt J <p(a)da. (III.11) 


k =1 1=1 


In order to simplify the notation in (III. 11) as well as facilitate the theoretical 
development that follows, we find it useful to define a parametric “information state” 
denoted by z v (t‘,a ) £ M, t £ [0,1], a £ A. For any a £ A, t £ [0,1], and set 
of searcher trajectories {x r,,k (s), 0 < s < f}, k = 1,2,. ..,K, z n (t\a) represents the 
cumulative detection rate given those searcher trajectories and vector of parameters 


a and is given by 


rt K L 


( EE r k,t (x v ’ k (s), y\s\ a)) ds, (III. 12) 

k =1 1=1 


or equivalently by the differential equation 


i ,? (s;a) = ^ (x^s), t/(s; a)) V s £ [ 0 , t], 


(III.13) 


k =1 /=1 


with ^(0;a) = 0. Using this notation, for any r] £ H, (III.11) simplifies to 

f(v)— [ exp (—2^(1; a)) <f> (a) da. (III.14) 

J a£A 

It is useful to simplify the notation in (III. 14) even further. We begin by 
defining the notation £ = (£ 1 ,..., £ K , 0) T £ M nA+1 , where f is a vector of initial states 
for the searchers such that the first K elements correspond to the initial “physical” 
states contained in rj and the final element corresponds to the initial “information” 
state. We then define the function F : M nA+1 x M nA " +1 — y M such that for any 
| £ M nAr+1 and x £ M nA+1 , where x = (x_i, z) T , with x_i £ M . nK , and z £ M, 


F(lx)±e~’. 


(III.15) 


To complete our notational simplification, for any a £ A, we also define the function 
/(•;«) : H —)■ M by 

f(Vioi) = F (i,x v (t;a)) , (HI-16) 
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where x v (t; a) is an augmented state defined by combining the “physical” states with 
the “information” state as follows 


r v (+'i 

^(t;a)= | * | eP ff+1 , 

z v (t; a) 


(III.17) 


where x v (t ) = (x v ’ 1 (t ) T ,..., x v,K (t) T ) T . Using this notation, for any rj G H the objec¬ 
tive function in (III. 14) simplifies to 


/(h)- / f(rj;a)<j>(a)da. 
JaeA 


(III.18) 


We now complete the definitions of ( GTP ) and ( GTP C ), which were prelimi¬ 
narily stated in Section II.D. Using the spaces defined in Section III.A we let 


and 


(GTP) min f(rj), 

H° 


(GTP C ) min f(rj). 

?;GH c 


(III.19) 


(III.20) 


2. Optimality Conditions 

In this section, we state our assumptions and give optimality conditions for 
(GTP) and (GTP C ). We begin by deriving parameterized differential equations of 
the augmented dynamics in terms of the augmented state, x(t ; a), defined in (III.17). 
For t G [0,1] and a G A, we define 

/ nM + \\ \ 


h(x(t),u(t ); a) = 


h 1 (x 1 (t ), u 1 (t)) 


nnK +1 


h K (x K (t),u K ( t)) 

\ZliZL !/'(«; a))/ 

where u(t) = (u l (t ) T ,..., u K (t) T ) T . Hence, 


(III.21) 


x(t] a) = h(x(t),u(t ); a), t G [0,1], £(0; a) = £. 


(III.22) 
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We note that x v f] •) is the solution of (III.22) when the input is rj = (£,u), and 
the augmented initial conditions are given by £ = (£ r , 0) T . We next state a series 
of assumptions, beginning with those related to (/>(•) and y l (•; •). It should be noted 
that in these assumptions and throughout the dissertation components of vectors are 
indicated by subscripts on variables. 

Assumption III.l. We assume that </>(■) is four times continuously differentiable. □ 

Assumption III.2. We assume that 

(i) for all t G [0,1], y l (t;-) is four times continuously differentiable for all l = 
1, 2,..., L, and 

(ii) y\-,a) is Lebesgue integrable on [0,1], for all a G A. 

□ 

The compactness of A and Assumptions III.l and III.2, respectively, ensure that the 
partial derivatives up to fourth-order of </>(•) and y l (t] •), t G [0,1], are bounded. The 
assumptions that </>(■) and y l (t\ •), t G [0,1], are four times continuously differentiable 
and the consequence regarding the bounds on their partial derivatives up to fourth- 
order are sufficient to complete the proofs of convergence in this chapter based on the 
later decision to use Composite Simpson’s rule for our numerical spatial integration 
scheme; see Section III.B.3b. If we had decided to use another numerical spatial 
integration scheme, such as Trapezoidal rule, the assumptions on </>(•) and y l (t\ •), t G 
[0,1], could be relaxed to being twice continuously differentiable with the consequence 
that their partial derivatives up to second-order are bounded. The assumption that 
y l (-,a) is Lebesgue integrable on [0,1], for all a G A is sufficient to ensure that the 
composite function r k,l (x k (-),y l (--,a)) is Lebesgue integrable on [0,1], for all a G A, 
under the assumption that r fc,/ (-, •) is continuous, and mild assumptions on x k f) as 
given in Assumption III.3; see Azagra et al. (2009). 

The next set of assumptions are related to r k,l (-,-) and /i(-,-;-). In these 
assumptions, and throughout the dissertation, we let h k (•, •) denote the n x n matrix 

Qi,fc / \ . 

of partial derivatives with element i,j given by , and let h k (-, •) denote the nxm 



Qljfc / \ 

matrix of partial derivatives with element (i,j) given by ^ ’ ; . We also adopt the 
matrix norm ||A|| = maxiuii =1 ||At;||, for any matrix A G M mxn , where v G M n is a 
column vector. We note that for matrices A G M mxn , B G M nxr , and column vector 
x G M n , we have that ||Ax|| < ||A||||x||, \\AB\\ < ||A||||.B||, and ||a;af r || = ||a;|| 2 (see for 
example p. 26 in Gill et ah, 1991). We also adopt the notation 

( h l( x { t )i u i t )) T \ 


h x (x(t),u(t)-,a ) = 


VEf=iE"=i V x r k ’ l (x k (t),y l (t,a)) T J 


h*(x(t),u{t)Y 

'L v-7 ^,k,l ( ry,k (+\ „,li 


(III. 23) 


where h x (x(t),u(t)] a) is a (nK + 1) x n matrix and 

( K{ x( A)w(t)) T \ 


h u (x(t),u(t)-, a) = 


hK(x(t),u{t)) T 

V 0 


( 111 . 24 :) 


where h u (x(t),u(t)] a) is a (nK + 1) x m matrix. 

Assumption III. 3. We assume for all k = 1,2,..., K and l = 1, 2,..., L that 

(i) there exists a C r < oo such that for all x k G M n and y l G M ny 


0 <r k \x k ,y l )<C r , 


(III. 25) 


(ii) r k,l (-,y l ) is continuously differentiable for all y l G M nj/ , 

(iii) r k ' l (x k , •) is four times continuously differentiable for all x k G M n , 

(iv) V x r k,l (x k , •) is four times continuously differentiable for all x k G M n , 

(v) h k (-,-) is continuously differentiable 

(vi) there exist C r i < oo, C r 2 < oo, and C r3 < oo such that for all x k G M n and 
y l G M ny , 

< C rl Vi = 1,2, j = 1,2,3,4, (III.26) 

< C r2 Vi = l,2,..., W , j = l,2,3,4, (III.27) 


d x r k,l (x k , y l ) 
dy\ 

di +1 r k’i(x k , y l ) 
dxdy{ 
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and 

\\V x r k ’ l (x k ,y l )\\ <C r3 , (III.28) 

(vii) there exists a constant K G [l,oo) such that for all x', x" G M n , v', v" G 
B(0, Pmax), and a G A, the following hold: 


|| h(x', v'; a) 

1 

IA 

>7 

7T 

1 

f~y 

+ 

7T 

1 

(III.29) 

|| h x (x',v';a) 

1 

H 

IA 

>7 

7T 

1 

+ 

tt 

1 

(III.30) 

and 



|| h u (x',v';a) 

1 

e 

IA 

>7 

7T 

1 

+ 

AT 

1 

< 5 ^ 

(III.31) 



□ 


We note that (III.29) implies that there exists a K' < oo such that for all x 1 G M n , 
and v G B(Q,p max ), 

||%»|| <lC[||x'|| + l]. (III.32) 

The assumptions about the detection rate function, r k,l (-,-), in Assumption 
III.3 are not overly restrictive as they allow for the use of many types of sensor 
models, and are similar to those used by other researchers (see Chung et al., 2010, 
for example). Assumptions III.3 (v) and (vii) are standard assumptions that parallel 
those adopted in Assumption 5.6.2 of Polak (1997). Assumption III.3 (vii) guarantees 
a unique solution to the differential equations governing the searcher dynamics given 
by (III.22). 

We next show that /(•) is Gateaux differentiable on H°, but in order to do 
this we need the following intermediate result. 

Lemma III.4. Suppose that Assumptions III.2 and III.3 are satisfied, then for any 
a G A, r] G H°, and Srj G 

(a) /(•; a) has a Gateaux differential Df(r /; cc; 5r)) at rj and hence a directional deriva¬ 
tive df(rj ; cc; Srj), with df(rj ; a; Srj) = Df(rj ; a; Srj), given by 

Df (ji\ or, Srj) = n f jrj; a), Sp) , (HI-33) 

' / H2 
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where the gradient 'V ri f(r}] a) = (V^/(t 7 ; a), 'V u f(r /; a)) T G -^ 00,2 is defined by 

V^/faa) = p"(0;o), (III.34) 

V u f(ir,a)(t) = h u (x v (t),u(t)]a) T p v (t]a), forte [0,1], (III.35) 


and jZ/t; a) is the solution of the adjoint equation 

P v (t] a) = -h x (x v (t),u(t)-,a) T p v (t-,a), t G [0,1], 

( 0 \ 


p\ l;a) = 


0 


g M nir+1 , 


(III.36) 


\ — exp (—^( 1 ; a)) J 


(b) the gradient V?,/(•; of) is Lipschitz H-continuous on bounded subsets of H° for 
all a E A with Lipschitz constant 

0 


K 


exp ( V2Ke K ) + K (V2Ke* + 1 


K 


(III.37) 


where K G [1, cx)) is as in Assumption III.3(vii) and K' < 00 is as in (III.32). 

(c) /(•;«) is Lipschitz H-continuous on bounded subsets of H° for all a G A with 
Lipschitz constant y/2Ke K , where K G [1, 00 ) is as in Assumption III.3(vii). 


(d) /(•; a) has a Frechet differential at 77 relative to H that is equal to . 0 /( 77 ; a 'i ^v)- 


Proof. Parts (a) and (d) result from an application of Corollary 5.6.9 in Polak (1997). 
Part (b) results from an application of Corollary 5.6.9 and the proof of Lemma 5.6.7 
in Polak (1997). For part (c), we deduce from Lemma 5.6.7 in Polak (1997) that 
z n { 1; a) is Lipschitz H-continuous as a function of ?/ with Lipschitz constant \/2Ke h . 
Since z v ( 1; a) > 0 for all 77 G H° and a G A, it follows that /(•;«) is Lipschitz 
H-continuous with Lipschitz constant \[2Ke K because the magnitude of the slope of 
the exponential function with an argument in the domain (— 00 , 0 ] is bounded by one. 
□ 

Proposition III.5. Suppose that Assumptions III.l, III. 2 and III. 3 are satisfied. 
Then, for any 77 G H (l and 5p G L 7 oo, 2 , /(•) has a Gateaux differential . 0 /( 77 ; a t V 

given by 

Df(trM = WMMn, ( IIL38 ) 

where the gradient V/(? 7 ) is given by 

Vf(rj)(t) = I \/J( m a)(t)<t>(a)dayt G [0,1], (IIL39) 

Ja&A 
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Proof. Let 5r) E H 00 , 2 , V £ H°, and a E Abe arbitrary. There exists a A* > 0 such 
that 77 + \5r] E H for all A E [0, A*]. For A E [0, A*], consider the ratio 


R(a ; 77 , 5rj, A) = 


/(t? + AJ 77 ; a) - f ( 77 ; a) 
A 


where from Lemma 111.4(a) we know that 


lim R(a; rj, Srj, A) = Df(ry, a; Srj) = (vf(ij-, a), 8r^j^ . 


From Lemma 111.4(c) we know that 


(III.40) 


(III.41) 


\R(a;rj,S V ,X)\ <L\\Sr}\\ H2 , (III.42) 

with L = \/2Ke K . We also know that 

V2Ke k [ (j){a)da = V2Ke k . (III.43) 

J aeA 

Then because Srj is a constant with respect to a, it is clear from (III.42) and (III.43) 
that for all A, /?(•; 77 , 5r), A) is dominated by an integrablc function. Then the Lebesgue 
Dominated Convergence Theorem yields that 


lim / R(a;rj, Srj, X)(p(a)da = / lim R(a] rj, 5r), X)<p{a)da. 


^4° Ja£A 


gA 


We then use (III.44) and (III.41) to obtain that 

Df(rrM = 1 lm Hn + M m a)-f(, r a) 
A^O A 


= lim 

A+O 


*eA 


f (77 + AJ 77 ; a) - /(77; q) 
A 


4>{a)da 


= lim / J2(q;; rj, Srj, A )(p(a)da 


^4° J a& A 

/ (V„/(77;a!),<S 77 ) 4>(a)da. 

JaeA N ' h 2 


(III.44) 


(III.45) 


Finally, ( V, ) /(t 7 ; a), Srj) is linear in Srj because by Lemma 111.4(d) /(t?; a) is Frechet 


772 


differentiable relative to H. Thus, (III.39) follows directly from (III.45). 
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Proposition III.5 also holds under weaker assumptions on </>(•), but because we 
need Assumption III. 1 later we adopt it here as well. This issue regarding Assumption 
III. 1 also applies elsewhere in this chapter. Our next task is to show that V/(-) is 
Lipschitz H-continuous on bounded subsets of H°. 

Lemma III. 6 . Suppose that Assumptions III.l, III.2 and III.3 are satisfied, then the 
gradient V/(-) is Lipschitz H-continuous on bounded subsets of H°. 


Proof. For any rf , rf £ H°, and t £ [0,1], 


l|V/(? 7 / )(t) — V/(r 7 ")(t)|| = / \/J{p';a)(t)(j)(a)da- W n ffr]"]a)(t)<p(a)da 

JaeA Ja£A 

From Lemma 111.4(b) we know for any a E A, V ,,/(•; a ) is Lipschitz H-continuous on 


bounded subsets of H n with Lipschitz constant e K 
Then, 


' exp ( 'V2Ke + K (\/2 Ke k + 1 


V^/(V; cx)(t)(j)(a)da — / V v f(r)"’, a)(t)(j)(a)da 


exp (y2Ke 


Because 


V/fa)(f) = 

we know from (III.46) that 


v £ /(u) 

v„/0)W 


1)] w - ■ 

(III.46) 

, t £ [ 0 , 1 ], 

(III.47) 


IIW/W) - V £ /(,,")H 2 + l|V„/W)W - V„/(r,")(f)|| 2 

< e*' exp (y2A'e A ') + A f/2Ke k + l) ||j/ - i/'Hj,.. (III.48) 

Hence, 

l|V/(77 / ) - V/(r 7 /, )||^ 2 

= l|V € /(77 / ) - V^/(r 7 ,, )|| 2 + f 1 \\VJW)(t)-V u f(ri /l )(t)\\ 2 dt 

Jo 

< e k ' exp (V2Ke k ^j + K (^/2Ke h + 1 j ||r/ — r/'H^ 

+ J e k ' exp (^/2Ke^ K (^/2Ke k + lj ||t/ — 77 "11 dt 

< 2 e k ' exp (V2Ke k ^J + K (V2Ke k + l) 2 b'- f/'IlL • (HI.49) 
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Because K and K' are both non-negative, (III.49) implies 

IIV/(r/) - V/(r/') \\ H2 < V2e k ' [exp (x/2 Ke k ) + K (V2Ke k + l)] \\rf - V "\\ H2 , 

(III.50) 

which completes the proof. □ 

We adopt the approach of Polak (1997) (see section 4.2), and state our optimal¬ 
ity conditions in terms of zeros of optimality functions; see also Section III.B.3a. We 
specifically define the nonpositive optimality functions 9 : H° —y M and 6 C : H c — y M 
as 

«(<?) = -i||V/()))|l5, 1 , (III.51) 

and 

° c (v) 4 min (V/fa), rj - rj) H2 + \ U - rj \\ 2 H2 , (III.52) 

trl c Zj 

which define optimality conditions for ( GTP ) and ( GTP C ), respectively. 
Proposition III.7. Suppose that Assumptions III.l, III.2 and III.3 are satisfied. 

(a) 9(-) and 9 C (-) are H° -continuous functions. 

(b) If rj G H° is a local minimizer of (GTP), then 0(f)) = 0. 

(c) If fj G H c is a local minimizer of (GTP C ), then 0 c ( r fj) = 0. 

Proof. The proof follows the same arguments as those for the proof of Theorem 4.2.3 
in Polak (1997), with Proposition III.5 taking the place of Corollary 5.6.9 from Polak 
(1997) and Lemma III.6 taking the place of Theorem 4.1.3 from Polak (1997). □ 

3. Consistent Approximations 

As discussed at the outset of this chapter, in order to solve (GTP) and (GTP C ), 
some form of discretization is necessary. In this section we define the approximat¬ 
ing problems ( GTP V), (GTPfi), (GTP NM ( N fi, and (GTPfi M(N ^). We also present 
conditions that must be satisfied to ensure that global minimizers, local minimizers, 
and stationary points of the approximating problems converge to global minimizers, 
local minimizers, and stationary points, respectively, of the original problems. We 
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divide our development into two subsections. Both subsections develop consistent ap¬ 
proximations for the pairs (( GTP),9 ) and (( GTP C ),9 C ), but the first subsection is a 
“stepping stone” that only deals with time discretization while the second subsection 
considers time and space discretization. In both subsections, we adopt the notation 
from sections 4.3 and 5.6 of Polak (1997). 


a. Time-Discretized Problems 

Let N denote the positive integers, and let Af be an ordered infinite 
subset of N defined by 


Af = {2 J }°^ 1 . (III.53) 

To begin our development, we must first define an infinite set of finite-dimensional 
subspaces H N C 77oo,2, whose union is dense in H 0 0j2 . For IV e Af and j = 0,1,.., IV — 
1, let the functions ttnj : [0,1] —» M, be defined by 

' VNWte[j/N, (j + 1 )/N), if j<N-2, 

TTNj(t) = VN vte [i/N , {j + 1)/IV], if j — N — 1, (III.54) 

0, otherwise. 

We also dehne the subspace WvCL™ 2 [0,l],by 


L n ± U6iy0,l] 


N -1 


«(■) = ^2ujKNj(-),Uj e M m ,Vj = 0,1,.., N — 1 

3=0 


(III.55) 


and let 


H n = M n xL N C (III.56) 

The functions 7 tnj(-) form a basis for L^ and are defined such that the relation be¬ 
tween Hn and the Euclidean space of coefficients used for numerical computation 
is isometric. We use the function space Hn for proofs of consistency of approxima¬ 
tion, and the Euclidean space of coefficients to develop implementable algorithms. 
The definition of the Euclidean space of coefficients, Hn, and further discussion of 
implementable algorithms are provided in Chapter V. We also dehne the sets 

Hjv = H n H n , (III.57) 
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(III.58) 


— H° n H n 

and 

H c> at ^ H c n H N , (III.59) 

with H. H° and H c defined in (III.5), (III.7), and (III.10), respectively. Based on the 
definition of A f, it is clear that the subspaces H N possess a desirable nested structure. 
This means that for any given N, N' G Af such that N' > N, Hn C Hn'- 

We will make use of Proposition 4.3.1 from Polak (1997), which is 
included here without proof for the sake of completeness. Below, we use the notation 
—4 to indicate convergence of a subsequence defined by A/". 

Proposition III.8. (a) Let H)) denote the closure of the set H°. Then, H”, —4 Hj), 
as N —>■ oo ? and (b) H r jv —4 H c , as N —y oo, where set convergence is in the sense 
of Painleve-Kuratowski 2 ; as N —» oo along the subsequence defined by Af. □ 


We now consider the approximate solution of (11.22) by means of for¬ 
ward Euler’s method. For any r/ = (£,u) G Hn and N G Af, we set xff^O) = f k and 
for any j — 0,1,..., N — 1, and for all k = 1, 2,..., K 

xt (0 + 1)/JV) - < urn = ifc* (xt(j/N), u\HN)) . (III.60) 

Simultaneously, we approximately solve (III. 13) also by forward Euler’s 
method. For any rj = (£,u) G H N , a G A, and N G Af, we set z v N { 0; a) = 0, and for 
any j = 0,1, ...,N - 1, 

4 ((i + «) - 4 (i4;«) = xf(j/N),y l (j/N-a)y (III.61) 

k =1 1=1 


Using the discretized “information state” given by the recursion (III.61), 
we define the approximate objective functions J'n '■ H^r —> M for any ^ G Hjv and 
N E Af by 

fN(rj) — [ exp (—44 ct)) (f)(a)da. (III.62) 

J aeA 


2 See Definition 5.3.6 in Polak (1997). 
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Again, for the sake of notational simplification, for any a G A, we also define the 
functions /#(•; a) : —» R. by 

In( iT, a) =F (i, x v N (t; a)) , (HI-63) 

where F(-, •) is as defined in (III.15) and x v N (j/N]a ) is an augmented state defined 
by 

iWi/iV; o) = ( 3:J ' (j/iV) ] 6 K"* +1 , j = 0,1, N- 1, (III.64) 

bjo/jv; a) y 

where x v N (j/N ) = (x v /(j/N ) T ,..., x v ^ K (j/N) T ^j , j = 0,1,..., IV — 1. Hence, for any 
IV G A/", we define the following approximating problems 

(GTP n ) min f N (n), (III.65) 

ven% 

and 

(CrTP c N ) min /at( ? ?)- (HI.66) 

j7£H c ,jv 

We note that the problems (GTPn) and (GTP^) still have spatial integrals that have 
not been discretized. 

We want to show that /n(-), N £ A/", are Gateaux differentiable and 
that their gradients are Lipschitz continuous. In order to do this, we begin with an 
intermediate result about the sums and products of Lipschitz continuous functions. 
We then prove an additional intermediate result which shows that the /jv(-; •) are 
Gateaux differentiable, and that their gradients are Lipschitz continuous. 

Lemma III.9. Suppose S is a bounded subset of H n . If F ,G : S —>• H° are Lipschitz 
H-continuous functions on S, then F + G, cF for any c G M, and FG are also 
Lipschitz H -continuous functions on S. 

Proof. The proof follows the same arguments as the proof of Theorem 4.6.3(b) in 
Sohrab (2003). In Theorem 4.6.3(b), however, the functions F,G are defined on a 
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bounded interval of M and our functions F, G are defined on a bounded subset of 
H u . Let L f and Lq be Lipschitz constants for F and G, respectively, then for any 
V, rf e S, 


\\(F + G)( V ) - (F + G)(t /)\\ H2 < l\F( V )-F(rf)\\H a + \\G( V )^G(rf)\\ Ha 

< (L r + L G ) \\v-v'\\h„ (III.67) 

and 

II (°F)(v) - ( c F)W)\\h 2 = \c\\\F( V ) - F{t]')\\h 2 < \c\L F \\r} - r/|| h 2 - (III.68) 

Since S is a bounded subset and F and G are both Lipschitz H-continuous on S, 
there exists M > 0 such that \\F(r])\\ H2 < M and ||G(77)||^ 2 < M for all rj G S. Then 
for any rj, 7 / G S we have 

\\(FG)( V ) - (fG)(t/)\\h 2 < \\G(v)\\h 2 \\F(v)-FW)\\h 2 

+ \\FW)\\h 2 \\G( V )-GW)\\ H2 

< {ML, + ML g )\\rf - rf\\ Ha . (IH-69) 

□ 

Lemma III. 10. Suppose that Assumptions III. 2 and III.3 are satisfied, and that 
N G A f, then for any a G A, 

(a) the function /jv(-;a) : —* M, is Gateaux differentiable, and the gradient 

V/v(h; a ) = (v^/jv(? 7 ;q;), V u f N {rj;a G H N is given by 

V ? /v(G«) = Pl(0;a), (III.70) 

N—l 

V u f N {rp,a){t) = '^2~i‘ n N {j/N;a)n N j(t), t G [0,1], (III.71) 
j =0 

where 

7jv0V-W; «) = ^ h u (x v N {j/N),u{j/N);a) T p n N {(j + 1)/N;a), 

J — 0, N—l, (III-72) 
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and //(,(•; a) is determined by the adjoint equation 


PnU / N ; «) - Pat((j + l)/iV; a) 

= jjhx(x v N (j /N),u(j/N); a) T p v N ((j + l)/iV; of), 

3 = 0 , 1 , — 1 , 

( 0 

Pat( 1 ;«) = q 

\-exp(-^(l;o;)) 

(b) For any bounded subset S of H. there exists a Lipschitz constant L s < oo such 
that, for any IV > 1, j = 0,1, N, 77 , 77 ' G <5 fl /7/v, and aei, 

/jv(»7; 01 ) - hirf; a) < Ls\\r) - p'\\ H2 , (III.74) 

and 

V/atC 77 ; «) - V/ Jv (77 , ; a) < L 5 ||t 7 - 77 'H^. (III.75) 

h 2 

Proof. The proof of part (a) is the same as the proof of Theorem 5.6.20 in Polak 
(1997), so it is not repeated here. For the proof of part (b), we begin by proving 
(III.74). From Theorem 5.6.16 in Polak (1997) we deduce that ^(l;a) is Lipschitz 
continuous as a function of 77 . Since z n N { 1; a) > 0 for all 77 G H° N and a G A, /jv(-; a ) 
is Lipschitz Hjv-continuous because the magnitude of the slope of the exponential 
function with an argument in the domain (— 00 , 0 ] is bounded by one. 

To show that V/n (S a ) is Lipschitz H^-continuous on bounded subsets, 
we look at its component parts beginning with (III.70). To show that V^/at( 77 ;q;) is 
Lipschitz Hjv-continuous on bounded subsets we proceed with an induction argument. 
From (III.70) we see that a ) is equal to the value of the adjoint at time zero, 

p v N ( 0; a). In order to find p v N ( 0; a), we use the recursion found in (III.73). We know 
that pir (1; a) is Lipschitz H^-continuous as a function of 77 for the same reasons that 
/jv(-;q;) is Lipschitz Hyv-continuous as a function of 77 . From Assumption III.3(vii) 
and Theorem 5.6.16 in Polak (1997) we also know that h x ( x v N (t) ,u(t) ; a) is Lipschitz 
Hjv-continuous as a function of 77 for any t G [0,1], and that its Lipschitz constant, 
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L~ hx < oo, is independent of t and N. Suppose that p n N ((j + l)/iV; a) is Lipschitz H r 
continuous as a function of rj and that its Lipschitz constant, L p < oo, is independent 
of N. Because h x (x v N (t) , u (t) ; a) and p v N ((j + l)/iV;a) are both Lipschitz Hat- 
continuous for all q e S, and S' is a bounded subset of H, there exists M < oo 
such that || h x ( x v N (t) , u (t) ; a) || < M and \\p n N ((j + l)/iV; a)|| < M. Then, based on 
(III.73) and Lemma III.9, we have that 


PnU / n ;«) - P%U/ N > a ) 


< 

+ 

< 


pn(U + 1 )/ Ar ; a ) ~ pn((J + 1 )/ N 'i a ) 

^ h x (x v N (j/N),u(j/N); a) T p v N ((j + 1 )/N] a) 
h x {x%{j/N),u(j/N)-,a) T p%((j + 1 )/N; a) 

L p\\v ~ rf\\ h 2 + + ML P ) \[q - 77 '|| h 2 

( m(^ + i)\ 

L P 1 + —H- L ||h - rfW (III.76) 


By doing another step in the backward recursion, we find that 




M[-p- + 1 

< L P | 1 +- h -- ) \\v - r/\\ h 2 


< L p I 1 + 


N 


M ( _|_ 1 

L/p 


N 


N 


\\v-v'\ 


h 2 - 


(III.77) 


Let Km = M y-j^ + 1J. There then exists an N e N such that for all N > N, 
(l + ^ l ) N < 2e KM . Then, for any N > N, we have 

Pl((j - 1 )/N; a) - pi((j - 1 )/N; a) < 2L p e K ™\\r 1 - rf\\ H2 . (III.78) 


For values of N smaller than N, we define 

N-i 


Si = 1 + 


K 


M 


N-i 


, i = l,2,..,iV-l, 


(III.79) 
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and 


K$ = max{2, max 5*}. (III.80) 

Then, 




<K s L p e KM \\ V - V '\\ H2 . 


(III.81) 


Hence, it follows by induction that V^/at(-;q;) is Lipschitz H^-continuous, and that 
the Lipschitz constant is independent of N. 

Next, we consider V u /at(-; a) given in (III.71). From (III.54) we know 
that 7 Ta Tj(t) is either 0 or y/N. The ^ factor in the definition of y^(-) in (III.72) 
combined with the 7r/vj(f) value of 0 or y/N ensures that V u /Ar(-;a) is independent 
of N. From Assumption III.3(vii) and Theorem 5.6.16 in Polak (1997) we find that 
h u {x v N (t) , u (t ); a) is Lipschitz H-continuous as a function of 77 for any t G [0,1], 
and that the Lipschitz constant is independent of t and N. Based on the induction 
argument given in (III.76) through (III.81), we know that p 71 N (^-;a) is Lipschitz 
LLv-continuous as a function of 77 for all j = 0,1, ...,N — 2, and that the Lipschitz 
constant is independent of N. From (III.54) we see that the summation in (III.71) 
is zero for all terms except where t is between j/N and (j + 1)/AL For these non¬ 
zero terms, V u /jv(s «) is composed of products of Lipschitz HAr-continuous functions 
whose Lipschitz constants are independent of N, so by Lemma III.9, V u /jv(? 7 ; a) is 
Lipschitz Hjv-continuous, and its Lipschitz constant is independent of N. This means 
that both components of V/jv( 77 ; a) are Lipschitz H^-continuous and their Lipschtiz 
constants are independent of N and a, so the proof is complete. □ 

Proposition III.11. Suppose that Assumptions III.l, III. 2 and III. 3 are satisfied, 
and N e Af, then for any 77 e and Srj G H^, }n{-) has a Gateaux differential 

Df N (irSrj) = (V/jv(? i),5v)h 2 > where 

V/Ar(? 7 )(t) = f V r; /jv(r/; a)(t)(j>(a)da, Vt G [0,1]. (III.82) 

JaeA 

Proof. The proof follows the same arguments as the proof of Proposition III.5 with 
/(•;•) replaced by /#(*; •), so it is not repeated here. □ 
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Next, we show that V/jv(-) is Lipschitz H^r-continuous on bounded 
subsets of H'y- 

Lemma III. 12. Suppose that Assumptions III.l, III.2 and III.3 are satisfied, then 
the gradient V/jv(-) is Lipschitz H N -continuous on bounded subsets of H^. 

Proof. The proof follows the same arguments as the proof of Lemma III .6 with /(•; •) 
replaced by /'v(•;■), so it is not repeated here. □ 

As in Section III.B.2, we state our optimality conditions in terms of 
zeros of optimality functions. For any N G A/", we define nonpositive optimality 
functions 9n ■ —> M and 0 C N : H C)JV —> M by 

Ml) & -i||V/ K ( f ,)|| 2 Hl , (III. 83) 

and 

&n(v) = niin (V/at(?/), 77 ' - rj) H2 + \ ||rf - 77 ^, (III.84) 

which characterizes stationary points of ( GTPn) and ( GTPf ), respectively. 
Proposition III. 13. Suppose that Assumptions III.l, III.2 and III.3 are satisfied. 

(a) 9]\r(-) and 9 C N (-) are H'y-continuous functions. 

(b) If i) G is a local minimizer of ( GTP N ), then 6 ^( 77 ) = 0. 

(c) If f/ G H Cj jv is a local minimizer of (GTPfj), then 6% ( 77 ) = 0. 

Proof. The proof follows the same arguments as the proof of Proposition 1.1.6 in 
Polak (1997), with the norms and inner products replaced with their H -2 equivalents. 

□ 

We are now ready to develop proofs for consistency of approximation for 
the pairs ((GTPn),6n) in the sequence {((GTPjsr), 9 n)}ngjv- Because we deal with 
more than one type of problem and its corresponding approximation, it is simpler to 
define consistent approximations and epi-convergence using abstract problems. We 
adopt the notation of Section 3.3.1 in Polak (1997) and let B be a norrned linear 
space, with norm || • ||g, and let M = M U {— 00 , + 00 }. We then define the problem 

(P) min f(x), (III.85) 
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where / : B — > M is lower semicontinuous and X C B is a constraint set. Then let 
{Bn}n&m be a family of finite-dimensional subspaces of B such that Bn\ C Bn 2 , for 
all Ni < N 2 G Af, and UngatBn is dense in B. For all N G N, let Jn ■ Bn —> K be 
a lower semicontinuous function that approximates /(•) on Bn , and let AW C Bn be 
an approximation to X. We then define the family of approximating problems 

(Pat) min Jn(x), N E Af. (111.86) 

xGXn 

Finally, we let A c / denote the closure of X and we define the problem (P c i) by 

(Pd) min f(x), (III.87) 

x&X c i 

which may not be epi-graphically equivalent to the problem (P), but which we will 
assume is equivalent to (P) in the sense that it has the same local and global minimiz- 
ers. Before we define consistent approximations, we must first define what it means 
for a function to be an optimality function. As on page 398 of Polak (1997), we define 
an optimality function as follows: 

Definition III.l. Let 5 be a subset of B such that X C S. We will say that a 
function 9 : S —» M is an optimality function for (P) if 

(i) 9(-) is sequentially upper semi-continuous, 

(ii) 9(x) <0 for all x <G S, and 

(hi) if a: is a local minimize!' of (P), then 9(x) = 0. 

Similarly, for all N G Af , let Sn be a subset of Bn D S such that Xn C Sn C S. We 
will say that a function 6n ■ Sn —>• M is an optimality function for (P/v) if 

(i) 6 n (-) is sequentially upper semi-continuous, 

(ii) 9 N (x) < 0 for all x G Sn, and 

(iii) if x n is a local minimizer of (P/v), then 6n(xn) = 0. □ 

Then, as on page 399 of Polak (1997), we define consistent approximations as follows: 
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Definition III. 2. Consider the problems (P), (P/v), and (Pd) defined in (III.85), 
(III.86), and (III.87), respectively. Let S' be a subset of B such that X C S and, in 
addition, Xn C S for all N G A f. Next, for all N G A7, let Sn be a subset of S' f) Bn 
such that Xn C Sn C S. Finally, let 0 a : S —» M and 9% : Sn —» M, A G A/”, be 
optimality functions for (P) and (P/v), respectively. We say that the pairs ((P/v), 0^), 
in the sequence {(Pn , 9%)} NeAT are consistent approximations for the pair ((P),9 a ) if 

(i) either (P/v) epi-converges to (P) or (P/v) epi-converges to (Pd), as A —» oo, and 

(ii) for any infinite sequence Ixa/Iats/c, AT C A/”, with a :at € S/v for all A G K, such 

that x n —>■ x as N —> oo, limsup^^.^ 9 < ^(xn) < 9 a (x). □ 

We base our proofs of epi-convergence on satisfying the conditions of 
Theorem 3.3.2 in Polak (1997), which we state here without proof for the sake of 
completeness. 

Proposition III.14. The epigraphs E N , N G A f, of the problems (P/v) converge to 
the epigraph E of the problem (P) if and only if 

(a) for every x G X, there exists a sequence Ix/vl/vgA/ - , with Xn G X n , such that 
x n -> A/ ” x, as A —* oo, and lirnsup Jn(xn) < /(x); and 

(b) for every infinite sequence {x/v}weA', with K C A7, such that x/v G AT^, for all 

A^ G K, and x/v —s, as A -> oo, i G I and liminf /jv(xjv) > /(x). □ 

We would like to establish epi-convergence of ( GTPn) to (GPP); un¬ 
fortunately, this is not possible. For this reason, we make a small modification to 
(■ GTP) and replace H° with its closure, H^. We define this new problem as follows 

( GTP d ) min f( v ). (III. 88 ) 

ven ° cl 

It is possible to establish epi-convergence of (GTPn) to (GTPd). We use the problem 
(GTP d ), with the following assumption. 

Assumption III.15. We assume that all local and global minimizers of (GTPd) are 

in H°. □ 

Using an approach similar to that found in Section 3.3 of Polak (1997), 
we next show that the pairs ((GTP N ),9 N ) in the sequence {((GTP n ),9n)}n£AT are 
consistent approximations for the pair ((GTPd), 9), which ensures that globally and 
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locally optimal points as well as stationary points of ( GTP N ) converge to correspond¬ 
ing points of ( GTP c i ), as N — * oo. In order to show the epi-convergence of ( GTPjy ) 
to (GTP c i), we will need the following intermediate result. 

Proposition III.16. Suppose that Assumptions III.l, III. 2 and III. 3 are satisfied. 
Then, given a bounded subset S C H, there exists a Ks < oo such that for every 
p G S fl Htv, arid N G A/", 

\Mv) - m\ < (III.89) 

Proof. For any p G S fl and N G A f, we have that 


\fN(v)~f(v)\ = 
< 


I aeA 


exp(—z^(l]a))(f)(a)da— / exp (— z v (l] a)) cj)(a)da 


'aeA 


'aeA 


|exp (— 2 ^( 1 ; a)) — exp (—z 1, ( 1; a))| f>(a)da. 


(III.90) 


By Theorem 5.6.23 in Polak (1997), we deduce that we can bound the approximation 
error (^(l; a) — z v N ( 1; a)|. Specifically, for any a G A, given a bounded subset ScH, 
for every p G S fl H/v, and N G A f, we have 

|«’'(l i a)-4(l;a)|<^, (III.91) 

with 

K s = max j 1, 1 + Kf, (III.92) 

and 

K^K'exp(k') supine'll+ 1 | (£,u)eS}, (III.93) 


K as in Assumption III.3(vii), and K' as in (III.32). Hence, by the properties of the 
exponential function and the fact that z]f(l;a) > 0 for all p G H u and z v (T,a) > 0 
for all p G Hjy-, for any a G A, and every p G S fl Hjv, and N G A f, 

l ex P (—2^(1; a)) — exp (—z v (l; a))| < |^(l;o;) - z v N (l;a)\ < (III.94) 

Because Ks does not depend on a, we have that 

I In(v) ~ fiv) I < / |exp 1; a)) - exp (-^(1; a))\<j>(a)da < (III.95) 

Ja&A 7V 

□ 


We now show that (GTPaj) epi-converges to ( GTP c i ). 
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Theorem III. 17. Suppose that Assumptions III.l, III.2, III. 3 and III. 15 are satis¬ 
fied. Then ( GTPn ) epi-converges to (GTP c i), as N —* oo. 

Proof. We show that ( GTPn ) epi-converges to ( GTP c i ), as N —» oo, by showing that 
the conditions in Proposition III. 14 are satisfied. Let rj € be arbitrary. Then, 
from Proposition III.8, there exists a sequence {r)N} N eAf such that rjN £ H^, for all 
N G A f, and rjN —^ r h as N —> oo. Let e > 0. By the H-continuity of /(•), there 
exists an TV G A/” such that for all N > TV, ^ < | and |/(^jv) — /('?/) | < where TVs 
is as in (III.92). Hence, by Proposition III. 16 

I/jvfav) _ /Ml < l/iv(w) - /(Vtv)I + |/Mv) - /Ml 

Ks e 
< — + - 
~ N 2 

< e, (III.96) 

for all N > N, N E Al. Consequently, /vC^/v) —M /M as At —» cx), satisfying 
condition (a) of Proposition III. 14. 

In order to show that condition (b) of Proposition III. 14 is satisfied, 
suppose that a sequence {t/tv }NeM su °h that 77 jv G for all N G A/”, and t/jv —M rj, 
as AT —>• cx). Then, based on the construction of in (III.58), we must have that 

rj G H/. It again follows from the H-continuity of /(•) and Proposition III.16 that 

/.v( r /v) —M / (77) as At —>• cx), which satisfies condition (b) of Proposition III.14. This 
proves that (GTPn) epi-converges to ( GTP c i ). □ 

In order to show that the pairs ((GTPn), On), in the sequence 
{((GTP N ),0 N )} NeA r, and the pairs (( GTP £), 0£r), in the sequence {((GTP^),0 c N )} N eM 
are consistent approximations for the pairs (( GTP d ),0 ) and (( GTP C ),6 C ), respec¬ 
tively, we will need the following two results. 

Proposition III. 18. Suppose that Assumptions III.l, III.2 and III.3 are satisfied, 
then for every bounded subset S C H° ; there exists a constant Kf < 00 such that, for 
any N G A/" and r) G S D Hn 

l|VM 7 ))-v/(i,)|| ffl <^. (111.97) 
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Proof. Based on the definitions of f(rj) and /at(?/) given in (III.14) and (III.62), 
respectively, and Propositions III.5 and III.11, 


' aeA 


V n f N (r]] a)(/>(a)da — / V v f(v, ot)<f>(a)da 


' a£A 


H 2 


< 


I aeA 


V v f N (v, «) - V„/(r r, a) <f>(a)da. 


\h 2 


(III.98) 


We now focus on 


V„/*fa;a) — V j) /(t 7 ;q;) . By Theorem 5.6.26 in Polak (1997), 


I h 2 


for any a G A, and for every bounded subset S C H°. there exists a Kp < oo such 
that, for all N G N, and p G S' D T/)v, 

K f 


V f/ /jv('7;«) - V„/(r7;o) 


< 

h 2 N 


(III.99) 


We deduce from the proof of Theorem 5.6.26 in Polak (1997) and Assumption III.3 
that Kp is independent of a. Then, we have that 



Vr Jn(v, «) 


V,,/(?/; a) 


4>(a)da < 

H 2 



(III.100) 


which completes the proof. □ 

Because we can only establish epi-convergence of (GTP/v) to ( GTP c i ), and not (GTPjy) 
to ( GTP ), we will need the following intermediate result. 

Lemma III. 19. Suppose denotes the closure of H". If a function V/3 : H° —> H 
is Lipschitz H-continuous on bounded subsets of H° ; then it is also Lipschitz H- 
continuous on bounded subsets of . 


Proof. Because H" is constructed from H by choosing 7 G (0,1) (see (III.7) and 
(III. 8 )), it is always possible to construct an H (,/ that is slightly larger than H 11 by 
choosing a 7 closer to one. Since H 1 contains H(), and it can be shown that 'Vfd(-) 
is Lipschitz H-continuous on bounded subsets of H lj/ , we can conclude that V/3(-) is 
Lipschitz H-continuous on bounded subsets of H|). □ 

We now show that the pairs ((GTP/v), 9n) in the sequence 
{((GTP/v), 0]\r)}NeAf are consistent approximations for the pair (( GTP c i),9 ). 
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Theorem III. 20. Suppose that Assumptions III.l, III.2, III.3, and III.15 are sat¬ 
isfied, and that ( GTP C \), 9, (GTPn), and On are defined as in (III.88), (III.51), 
(III.86), and (III.83), respectively. Then, the pairs (( GTPn ), On), in the sequence 
{(( GTP n ), On)}ngjv are consistent approximations for the pair (( GTP c i),9 ). 

Proof. Suppose that an infinite sequence {'On} NeN suc ^ that t)n £ Hjy, for all 
TV e A f, and t]n — ^ V, as TV — y oo. Let e > 0. From Lemma III .6 and Lemma III. 19 
we know that V/(-) is Lipschitz H-continuous on bounded subsets of H° ; . By the 
H-continuity of V/(-), there exists an TV e N such that for all TV > TV, < | and 
II V/(r? 7 v) — Vy( 77 )|| < |, where Kp is as in Proposition III.18. Hence, by Proposition 
III.18 

I|V/jv(i7jv) - V/(i7)||/f 2 < ||V/ iV M-V/M||^ + ||V/M-V/(r 7 )||H 2 

K F e 
< — + - 
“ TV 2 

< e, (III.101) 

for all TV > TV, TV e J\f . Consequently, V/jv(77jv) —V/( 77 ), as TV —y 00 , implying 
that On(t]n) — 0(rj), as TV —> 00 . Theorem III.17, together with the convergence of 
On {'On) to 6(rf) as TV —» 00 , satisfies the requirements of Definition III.2 for consistency 
of approximation. □ 

Next, we show that the pairs ((GTPff) j On) in the sequence 
{((GTfft), are consistent approximations for the pair (( GTP C ),6 C ). 

Theorem III. 21. Suppose that Assumptions III.l, III.2, III.3, and III. 15 are sat¬ 
isfied, and that ( GTP C ), 9 C , (GTPff), and 0 C N are defined as in (III.20), (III.52), 
(III.66), and (III. 84 ), respectively. Then, the pairs ((GTPff), 0 C N ), in the sequence 
{{(GTP c n ), O c n )}n&m, are consistent approximations for the pair (( GTP C ),9 C ). 

Proof. The proof that the problems (GTPff) epi-converge to ( GTP C ) is the same as 
the proof of epi-convergence of the problems (GTPn) to (GTP c i) given in Theorem 
III. 17 above. 

From the proof of Theorem III.20 above, we know that V/at(t 1 /v) — 

V f(rj), as TV —> 00 . Then, following the same arguments as in the proof of Theorem 
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such that 


4.3.6 in Polak (1997) we see that given any infinite sequence {vn}n & j^, 
r/N £ Hc.n for all N £ Af, which converges to an rj £ H c , 9 c n {j]n ) — 9 c (i 7 ), as IV —* 
00 . Theorem III.20, together with the convergence of 9 c n (t]n) to 9 c (r /) as iV —* 00 , 
satisfies the requirements of Definition III.2 for consistency of approximation. □ 

b. Time- and Space-Discretized Problems 

We next consider the time- and space-discretized problem. As discussed 
in Section II.A we focus on A C M 2 , and introduce the space discretization parameter, 
M = (Mi, M 2 ) T £ N x N. We also define a generic numerical integration scheme, Im, 
for a function, T : C P (A ) —» M, where C p represents the differentiability class . 3 The 
integration scheme Im is a mapping from C P (A ) to M. or to M ri+m for any M eNxff 
and T £ C P (A ) defined by 

Mi M2 

j=l j =1 

where i = 1 , 2 , ,.,, Mi, j = 1 , 2 , ...,M 2 , are weights for the chosen numerical 
integration scheme, and a tJ are the discretization points at which the integrand is 
evaluated. 

We then make use of the integration rule, I M , to define the approximate 
objective function f^M ■ —y M for any rj £ H^, N £ A I, and M 6 N x N by 

fNAfiv) ~ (exp [-z v N ( 1; •)] </>(•)) • (III.103) 

We next consider the differentiability of /atm( - )- 

Proposition III. 22. Suppose that Assumptions III.l, III.2, and III.3 are satisfied, 
Im is defined as in (III. 102), N £ Af, M £ N x N, and f^M '■ H^r —> M is defined as 
in (III. 103), then for any rj £ H° N and Sr/ £ fNAi(-) has a Gateaux differential 

Df NM {wM = (V fNM(ri),$o) H2 , where 

Vf N M(v)(t) = Im [vJ N { m •) (t)<f>(-)] , Vt £ [0,1], (III.104) 

3 The class C° consists of all continuous functions on A. For any positive integer, p, C p is the set 
of all differentiable functions on A whose gradient is in C p_1 . 


T(q;) da , 


(III.102) 
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Proof. Let Sr) E and r) E H^r be arbitrary. From Lemma 111.10(a) we know 


Df N (rr, or, Sr)) = ( V v f N (v’, «)> 


h 2 


(III.105) 


Then, we have 


DfNM^Srj) = 


f nm ij) + AfSry; a) - f NM (m a) 


lim 

A4.0 

Mi m 2 

‘IjsEE 1 ^ 

1=1 j =i 

Mi M 2 

E E w « “s 

*=i j=i 

Mi M 2 


A 


/tv (?7 + XSr)-, on ,) - /at(77; « 



E \^vfN(V', oiij), Sr)) 0(aij) 

11 ' ' H2 

1=1 j=\ 

I Mi M 2 \ 

EE W « V 

*=1 t =1 / H 2 


(III.106) 


Our next result is related to the Lipschitz Hjv-continuity of V/atm(^) 
on bounded subsets of H^. 

Lemma III.23. Suppose that Assumptions III.l, III.2, arid III.3 are satisfied, then 
the gradient V/jvm(^) A Lipschitz H ^-continuous on bounded subsets ofH%. 


Proof. The proof follows the same arguments as the proof of Lemma III. 12, with 
integration replaced by Im- □ 

In the analysis that follows, it is necessary to quantify the error in¬ 
troduced by the numerical integration scheme in order to complete the convergence 
proofs. Clearly, the choice of numerical integration scheme determines the relationship 
between the discretization level and the amount of error introduced by the approxi¬ 
mation. In order to conduct our analysis, we make the following assumptions. 
Assumption III.24. We assume that there exist scalars a,b,c, and d such that 

(i) A = {(a 1 ,a 2 ) £ M 2 |a < aq < b, c < a 2 < d}, where A represents the region of 
integration in Im- 
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(ii) aij = (ai(i), a 2 (j)) T and an integer Mi is chosen so that the interval [a, b] is 
subdivided into 2M\ subintervals {[aqty — 1),Q'i(i)]}^ 1 of equal width h = 
(i b — a)/2Mi by using the equally spaced discretization points aq(i) = a + ih 
for % = 0,1,2, 2M\. An integer M 2 is chosen so that the interval [c,d] is 
subdivided into 2M 2 subintervals {[a 2 (j — 1 ), « 2 (j)]}j=i 2 of equal width k = 
(d — c)/2M 2 by using the equally spaced discretization points a 2 (j) = c + jh for 
j = 0, 1 , 2 , ..., 2 M 2 . Composite Simpson’s rule is used for numerical integration 
and for T G C P (A) is given by 


T , hk 

hi{^) - — 


m 2 


d/(a, c) + T(a, d ) + T( b,c ) + T(5, d) + 4^ T(a, a 2 (2j - 1 )) 

j =i 

M2 — 1 M2 M2 — 1 

2 ^ ^(a,a 2 (2j)) + 4^^(6,a 2 (2j-l)) + 2 d/(6,a 2 (2j)) 

l=i l=i l=i 

Mi Mi-1 Mi 

4 d / (o;i(2?'— 1), c) + 2 T(aq(2i), c) + 4 T(ai(2i — 1), <f) 

i= 1 2=1 2 — 1 

Mp —1 M 2 Mi 

2 ]T ^(a 1 ( 2 i),d) + 16^ ^^( ai ( 2 i-l),a 2 ( 2 i-l)) 

2—1 J = 1 2—1 

M 2 -l Mi M 2 Mi-1 

*ee d , (a 1 (2i — 1), a 2 (2j)) + 8 EE d/(o;i(2i), a 2 (2j - 1)) 

3 =1 *=i j'=i *=i 

M 2 - 1 M 1-1 \ 

*E e v h(a 1 (2i), a 2 (2j)) j (III.107) 

i=i i=i / 


□ 


We note that the convergence proofs given below in Proposition III.26 would follow 
similar arguments if A had higher dimensionality, or a different numerical integration 
scheme had been utilized. The proofs could also be done if A was assumed to be a 
shape other than rectangular, but they would be more complicated. 

We find it necessary to show that the partial derivatives of faiv ; •)(/>(') 
and V f/ /jv( 77 ; •)</>(•) up to and including the fourth-order are bounded for any rj G 
a G A, and N G J\T in order to complete the proofs of convergence of /jvm(»)) to 
f{rj) and V/vm(^) to V/(r/), based on the choice of Composite Simpson’s rule as the 
numerical integration scheme. To facilitate these proofs, we begin by defining some 
notation. For any rj G H 1 ^, a G A, N G A/”, and j = 0,1,..., N — 1, we define 
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C 3 (a)=^(0;a), (III.110) 

and 

Ci(«) = Pn b • (III.lll) 

We note that Ci(')> Cs(')> an d C 4 (') depend on p and N. We next show that the partial 
derivatives of Ci(')> •••,C 4 (') up to and including the fourth order are continuous and 
bounded for any p G H^-, a G A, and N G A/". 

Lemma III. 25. Suppose that Assumptions III.l, III.2, and III.3 are satisfied and S 
is a bounded subset ofH () N . Then, 

(i) G(-)eC 4 (A),* = 1,2,3,4, 

and 

(ii) there exists C < oo, such that for all p G S, j = 0,1, ...,N — 1, a G A, and 
N G M 

<C Vi = 1, 2, V/i = 1, 2, 3,4, Vk = 1, 2, 3,4. (III.112) 

Proof. It can be seen by repeated applications of the chain and product rules that 
based on Assumptions III.2(i) and III.3(iii), Ci(') G C 4 (A). We also know from As¬ 
sumption III.l that C 2 ( - ) ^ C 4 (A). 

We now consider ( 3 (a). We know that p’ 1 N ( l;a) G C 4 (A) for the same 
reasons that Ci(-) G C 4 (A), because p n N (l] a) is a column vector of nK zeros and — Ci( - )- 
We also know that h x (x v N (j,), u(^); a), j = 0,1, ...,N-1, has h k (x v N {jj),u{jj)) , k = 
1, 2, ...,K, and Ef=i Et 1 ( x n( n)i y l ( n'i °0) as components. Because there is 

no a dependence, h k and partial derivatives of h k (x v N (fi), with 


da-i 
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respect to a are constants for all k — 1,2, Based on Assumptions III.2(i) and 

III.3(iv), V x r k ’ 1 •)) G C 4 (A). This means that h x (^(^), «(^);') e 

C 4 (A). Next, we proceed with an induction argument. Suppose that p^((j + l)/iV; •) G 
C 4 (A), then p v N (j/N; •), given by 

PnU/ n > •) = Pvr((i+ 1 )/ Ar ; + n(j/iV); •) T P^r((i+ 1 )/^; ■), (HI.113) 


is also C 4 (A). Hence, it then follows by induction that p r j^(0 ; •) G C 4 (A), and therefore 

0.(0 G C\A). 

By the induction argument above, we know that p v N ( bbl; •) G C 4 (A) for 
all j — 0,1,..., N — 2. Hence, £ 4 (-) G C 4 (A), and the proof of part (i) is complete. 

We start the proof of part (ii) by looking at Ci(°0- If can be seen by 
repeated use of the product and chain rules that the first through fourth-order partial 
derivatives of C^-) are made up of sums, differences, and products of the expressions 

»p [-x^£X> M (a (£) y (7“))1 (IIL114) 

. j= o k=i i=i v v 7 v 7 7 _ 

and 




Vi = 1,2 ,Vk = 1,2, 3,4. (III.115) 


To prove that the first through fourth-order partial derivatives of (j(-) are bounded, 
and that the bounds are independent of N and rj, we show that expressions (III. 114) 
and (III. 115) are bounded with respect to a, and that their bounds are independent 
of N and rj. We first consider (III. 114). By Assumption III.3(i), there exists C r < oo 
such that 

r kt (*t (A) J ( 7 “)) | - C ” (IIL116) 
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for all r) G S, k = 1, 2, K , l = 1, 2,..., L, a G A, j = 0,1 ,N — 1, and N G A/”. 
Then 


JV—1 K L 


N 

j= 0 fc=l z=i 

iV-1 AT L 

* vEEE 

j=0 fc=l i=l 
, AT-1 


AT 


i / •? 


^ u? >2/ 


AT 




AT 


< — V KLC r = KLC r 
- N ^ 


(III.117) 


3=0 


Hence, (III. 114) is bounded by exp(— KLC r ), which is independent of N and rj. 

Next, we consider (III.115). By Assumptions 111.2(h) and III.3(vi) 


d h 

da? 


k,l ( „v,k ( l_ 


r ’ \ X N 1 ft 1 rV l 


1 A 


AT 


< C r \Cy , 


(III.118) 


for all rj G S, k = 1, 2,..., K, l = 1, 2,..., L, a G A, i = 1, 2, j = 0,1 ,..., N — 1, 
k = 1,2,3,4, and N G Af, with C' r i and C y as in Assumptions 111.2(h) and III.3(vi), 
respectively. Then, for all i = 1,2, and n — 1, 2, 3,4 


N— 1 ^ K L 


A - <9af 

j=0 1 k= 1 i=l 


AT’ 


(III.119) 


< 


< 


N-l K L 

vEEE 

j= 0 k= 1 i=l 
iV-1 

- ^ KLC rl C y = KLC rl C y . 



(III.120) 


3=0 

Therefore, (III. 115) is bounded by KLC r iC y , which is independent of N and rj. Be¬ 
cause (III.114) and (III.115) are bounded, and their bounds are independent of N 
and rj, the first- through fourth-order partial derivatives of Ci(«) with respect to a 
are bounded, and the bounds are independent of N and rj. 

Now, we consider ( 2 ( 0 ;). From Assumption III. 1 and the compactness 
of A, we know that the first- through fourth-order partial derivatives of £ 2 (a;) are 
bounded for any a G A. 
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We next consider £ 3 (a). We know that the first- through fourth-order 
partial derivatives of p^( 1; a) with respect to a, are bounded for the same reasons 
that the first- through fourth-order partial derivatives of £ 1 ( 0 ;) are bounded. The 
first- through fourth-order partial derivatives of h k r (x^(-), «(•)), k = 1,2, with 

respect to a are bounded because they are all equal to zero. The first- through fourth- 
order partial derivatives of V x r k ’ 1 [x n N {t) 1 y l {t] a)) with respect to a are bounded for 
any t E [0,1] by Assumptions 111.2(h) and II1.3(vi), and the bounds are independent 
of t and N. Then, because h x (x(t),u(t ); a) is earlier defined by 

( K.(x(t),u(t)) T ^ 

h x (x(t),u(t);a) = „ , (111.121) 

hf(x(t),u(t)) T 

(E%=i Eti V x r k ’ l (x k (t),y l (t; a)) T j 
there exists K~ < oo such that 

< K' k ,, (IH.122) 

for all y E S, a e A, i — 1, 2, j — 0,1,..., AT — 1, k = 1,2, 3,4, and N E A f. 
From Assumptions III.3(ii), (iii), and (v), we know that h x (x v N {j/N),u(j/N);a) is 
continuous with respect to a for all rj E S, on A. Since h x (x v N (j/N), u(j/N)] a) is 
continuous with respect to a, and S is a bounded subset of H^, there exists < oo 
such that 

\h x (x^ N (j/N),u(j/N) ]a )\\ < K~ hx , (111.123) 

for all rj E S, a E A, j = 0,1,..., N — 1, and N E A f. Because zj^( 1; a) > 0 for all 
rj E S and a E A, there exists C\ < oo such that ||p^(l;o;)|| < C\ for all rj E S, 
a E A, and N E AT. To show that p v N (j/N ; a) is bounded for all y E S, a E A, and 
j = 0,1,..., N — 1, we proceed with an induction argument. Suppose that there exists 
Cj + 1 < oo, such that \\pjf((j + 1 )/lV; a)|| < C ]+ \. for all y E S, a E A, and N E A f. 
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Then, based on (III. 113) 


WPnU/N'i a) || < \K((j + l)/N;a) 

h x (x v N (j/N),u(j/Ny, a) T p v N ((j + 1 )/TV; a) 


TV 


< C j+ 1 + —K- hx C j+ 1 - C j+ 1 


1 + 


A r 


AT 


(III.124) 


By doing another step in the backward recursion, we find that 


I \Pn(U - l)/iV; «)|| < C j+1 1 + 


Kt 


TV 


+ 


Kt 


N 


-C 


3 +1 


1 + 


Kt 


N 


= C 


3 +1 


- Cj'+l 


1 + 


K h ' 2 

'I'X 

~w 


1 + 


Kt 


N 


N 


(III.125) 


There exists IV e N such that for all N > TV, ^1 + < 2e Kh *. Then, for any 

TV > TV, 

IK((i - 1)/N; a)|| < 2C i+1 e^. (III.126) 

For values of TV smaller than TV, we define 

/ k~ \ 

+ , i = l,2,...,TV-l, (III.127) 

and 

Kg — max{ max {<5j},2}. (III.128) 

i=l,2,...,AT-l 

Then, 

WpKU - 1)/A r ;a)|| < K s C j+1 e K -.. (III.129) 


Hence, it follows by induction that p v N (j/N]a ) is bounded for all p S, a A, 
j = 0,1,..., N — 1, and the bound is independent of N and r\. 

We now proceed with another induction argument to show that the 
first- through fourth-order partial derivatives of p v N (j /TV; a) with respect to a are 
bounded. Suppose the first- through fourth-order partial derivatives of p n N ((j + 
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1)/N;a) with respect to a are bounded and that the bound, C') +1 < oo, is inde¬ 
pendent of N and r/. Then, based on (III.113), for all rj G S, a G A, i = 1,2, 
j = 0,1,..., IV — 1, k = 1, 2, 3,4, and N £ J\f 


d K 


1 

'N 


da K 

d K 


P V N (j/ N 'i a ) 


< 


d K 


:Pn(U + i-V-W; a) 


da? 


da** 

h x (x v N (j/N),u(j/N); a) T p v N ((j + 1 )/N; a] 


< h +1 + ^ 


1 

N 


d K 


h x (x v N (j/N),u(j/N); af—p'M + 1 )/N\ a 


hx(x v N ti/ N ),u(j/N); a) T p v N ((j + 1 )/N; a 


< ^ + i + 


- I 1 + 


<9cy 

K~ h C' j+l + C j+ 1 K'~ h 


N 


Kr 


+ 


K'tC j+1 


N J N 

By doing another step in the backward recursion, we find 

d 


(III.130) 


dal 


Pn(U - 1 )/ N ~>«) 


— Cj + 1 I 1 + 


— Cj + 1 f 1 + 


Kt V KLC j+1 _ C, +1 | C,A3 


AT 


+ 


AT 


+ 


N 2 


+ 


N 


K, \ N KL Cj+1 K~ h K>~C j+1 CjK'~ 


N 


+ 


N 


+ 


N 2 


+ 


< K s C' +l e K ^ + K c , 


N 

(III.131) 


where Kc < oo is such that 


K'~ C j+1 K h K- C j+1 CjK'r 

A C > -77-t-777-h 


N 


N 2 


N 


(III.132) 


for any N e A f. Hence, it follows by induction that the first- through fourth-order 
partial derivatives of p\ r (0; •) with respect to a are bounded, and therefore the first- 
through fourth-order partial derivatives of Q$ (•) with respect to a are bounded. Fur¬ 
thermore, these bounds are independent of N and rj. 

By the induction argument above, the first- through fourth-order par¬ 
tial derivatives of p n N ( •) with respect to a are bounded for all j — 0,1,..., N — 2 
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and the bounds are independent of N and rj. This means that the first- through 
fourth-order partial derivatives of £ 4 (-) with respect to a are bounded, and the proof 
of part (ii) is complete. □ 

In order to prove epi-convergence and consistency of approximation, we 
need the following proposition. 


Proposition III.26. Suppose that Assumptions III.l, III.2, III.3, and III.24 are 
satisfied. Then for every bounded subset S C H° ; there exist constants Kj\ , Kj 2 < oo 
such that, for any N e J\f, M e N 3 x N 3 , where N 3 = {m e 2N + l\m > 3}, and 

tj g s n H]\f , 


and 


I/(h) - f nm{j ])| < 


Ks K n Kn 

N + (Mi - l) 4 (M 2 - l) 4 


(III.133) 


||V/(h)-V/vM(h)lk< 


Kf K n K n 

N + (. M 1 - l) 4 (M 2 - l) 4 ’ 


(III.134) 


and K s and K F are defined as in (III. 92) and Proposition III. 18, respectively. 


Proof. We know that 


I/(h) - /iVAf(h) I < I/(h) - /iv(h)l + |/jv(h) - f nm(j]) I, (III.135) 


and 


IIV/(h)■ V/jvAf(h) 11^2 < l|V/(h)-V/ A r(h)||^ + ||V/ i v(h)-V/ i vM(h)IU 2 . (111.136) 

Based on Proposition III.16 and Proposition III.18, respectively, |/(h) — /jv(h) I — 
and || V/(h) — V./v(h)II, for any r) G 5 fl %. In order to deal with 
IMh) - /jVAf(h)l and II V/iv(h) - V/iVM(h)ll / 7 2 , we define for notational simplicity 

<7i(a) 4 Ci(a)C 2 (a), (III.137) 


and 


92 (a) = V/v(h;«)/(«), 


(III.138) 



where by Lemma III.10 and the definitions given in (III.109), (III.110), and (III.Ill) 
we can write out the components of ( 72 (a) as 

g 2i (a) = V ? /iv(r 1 ; a)<j>{a) = p v N ( 0; a)<j>(a) = C 3 (a)C 2 (a), (III.139) 


and for t e [ 0 , 1 ] 
92ut(oi) 


N—l 


= V u /jv(77;a)(t)0(a) = 

K {A Of) ,«(£)) 


( ^;a) n N j(t) 
3 =0 

T 

( hi (rrV (3_\ (L\\\ 


</>(«) 


JV-1 

E 

3=0 


h u W (&)>“(£)) 
0 


vr 7 Vj (t)C4(a)C2(a). (III.140) 


Then, by Lemma III.25(i), g\(-) E C A (A) and g 2 (-) £ C 4 (A). Again, for notational 
simplicity, we define 


E\ 4 
^ 2 * = 


Im(9i(-))~ / <71 (a) da 

J aeA 


ImAA)) ~ / g- 2 &(a)da 

J a£A 


(111.141) 

(111.142) 


and 


E 


A_ 

2 utj — 


lM{g2utj{-)) 



g 2 utj{ot)da , 


(III.143) 


where g 2 & (a) is the i th component of g 2 % (a) and g 2 utj(a) is the j th component of 
<? 2 ut(a). From Lemma III.25(ii) for p — 1, 2£i, 2utj, i = 1,2,..., n/l + 1, t e [0,1], and 
j = 1, 2, ..., m, there exist constants Ci,C 2 < 00 that have no dependence on 77, a , or 
At such that 


max 

aeA 


dV(a) 

<9af 


< 6 ', 


(III.144) 


and 


max 

aeA 


d 4 g P (a) 
da 2 


< C 2 . 


(III.145) 
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Then, under Assumption III.24, when using Composite Simpson’s rule to approximate 
the integral of gi(-), E\ is bounded by (see, for example, pages 127-128 in Faires & 
Burden, 1993) 



(b — a) 4 d 4 g 2 (Li {a) (d - c) 4 d 4 g 2 ^(a) 

(Mi - 1) °eA daf (M 2 - l) 4 «eA da\ 

(d - c)(b - a) (b - a) 4 (d - c ) 4 < Kn K n 

180 [(Mi - l) 4 1 (M 2 - l) 4 ' 2 J “ (Mi - l) 4 + (M 2 - l) 4 ’ 


E 2 £i < 


(d — c)(b — a) 
180 


where Kn and Kj 2 are as above. Similarly, under Assumption III.24, when using 
Composite Simpson’s rule to approximate the integral of g 2u tj(-), E 2ut j is bounded for 
t G [0,1] and j = 1, 2,..., m by 


(b — a) 4 d 4 g 2ut {a) (d - c) 4 d 4 g 2ut (a) 

(Mi - l ) 4 «eA daf (M 2 - l ) 4 daf 

(d - c)(b - a) f ib-a) 4 (d - c) 4 1 < K n K r > 

180 - l) 4 1 + (M 2 - l) 4 2 J “ (Mi - l) 4 (M 2 - l) 4 ’ 


c, / (d-c)(b-a) 

2utj ~ 180 
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where, again, K n and Kj 2 are as above. Then, we have 


< 


IIV/wM - V/™(,) f H2 

l|V 5 / K (r;)-V { / NM (r;)|| 2 + [ \\VJ K (r,)(t)-VJ N „(r,)(t)\\ 2 dt 

Jo 

"v ( Kli V Kn Kn V 

h \( Mi - x ) 4 ( M2 - x ) 4 / h v( Mi - x ) 4 ( M2 - i) 4 / 


(nK + 1 + m) 


K 


n 


K 


12 


(Mi - l) 4 + (M 2 - 1) 


(III.148) 


We let K n = \fnK + 1 + mKn and Jl /2 = V'nK + 1 + rriKf 2 , and (III. 134) follows 
from (III. 136). Finally, 


AT, K T2 ^ Kn Ki 2 

~ (Mi - l ) 4 + (M 2 - l ) 4 “ (Mi - l ) 4 + (M 2 -l) 4 ’ 

(III.149) 

and (III. 133) follows from (III. 135), which completes the proof. □ 

For any N e A/", and M e N x N, we dehne the following approximating 

problems 


(GTPmm) 

min 


(III.150) 


(GTP c nm ) 

min 

V 6H c ,jv 

/jvm(j))' 

(III.151) 


For any N e A/", and M e N x N, we dehne nonpositive optimality 
functions (Tvm : H ( ( r —* M and d c NM : H Ci at —* M by 


^m(a('))- / gi(a)da 


I a£A 


Onm(v) ~ f NM iji)\\ 2 H^ 


(III.152) 


and 

^m(^) - nhn (Vf NM (v),v' - h)// 2 + J IIV - hll^ 2 , (III.153) 

T] G xl C) jV ^ 

which characterize stationary points of ( GTP^m) and (GTP^ M ), respectively. 

Proposition III.27. Suppose that Assumptions III.l, III.2, III.3, and III.24 are 
satisfied. 
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(a) 0 NM (-) and 0 C NM (-) are H^-continuous functions. 

(b) If fi G H” v is a local minimizer of ( GTPnm ), then Onm(v) — 0- 

(c) If i) G H Cj jv is a local minimizer of (GTP^ M ), then 0 c NM (f)) = 0. 

Proof. The proof follows the same arguments as the proof of Proposition 1.1.6 in 
Polak (1997), with the norms and inner products replaced by their H -2 equivalents. □ 
The proofs of convergence that follow require that we have only one 
discretization parameter, therefore for % = 1,2, we define M t : A/" —* N and make the 
following assumption. 

Assumption III.28. We assume for i — 1, 2 that M^(N ) — * oo, as N — > oo. □ 

It is not possible to establish epi-convergence of ( GTP NM ^ N )) to (GTP), 
but it is possible to establish epi-convergence of (GTP NM (n)) to ( GTP d ). We now 
show that the family {((GTPnm(n)), 0nm(n))}ngjv is a sequence of consistent ap¬ 
proximations for (( GTP c i),9) by showing that the conditions of Definition III.2 are 
satisfied. 

Theorem III.29. Suppose that Assumptions III.l, III.2, III.3, III. 15, and III.28 are 
satisfied, ( GTP d ), 6, (GTP NM ^ N )), and d NM ^ are defined as in (III.88), (III.51), 
(III. 150), and (III. 152), respectively. Then the pairs (( GTP NM ^ N )), 0 NM ( N ^), in 
the sequence {((GTP/vvqv)), 9nm(n))}n£AT are consistent approximations for the pair 
({GTPd), 0). 

Proof. We first show that (GTPvm(v)) epi-converges to (GTPd), as N —> oo, by 
showing that the conditions in Proposition III. 14 are satisfied. We begin by showing 
that part (a) of Proposition III. 14 is satisfied. Let p G be arbitrary. Then from 
Proposition III.8 there exists a sequence {Vn}n g j^ such that r/jv G H^, for all N G A f, 
and r/N —P v " p as N —> oo. Let e > 0. By Assumption III.28 and the H-continuity 
of /(•), there exists an N G Af such that for all N > N, i) 4 — f > 

^ < |, and \f(pN) ~ f(v) I < where Kj\ and K I2 are as in Proposition III.26, and 
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Kg is as in (III.92). Hence, by Proposition III.26(i), 


| fNM{N){m) ~ f(v )I < | fNM(N)(m) - f(m)\ + I/M - f(v )I 

Kn K I2 e 

~ N + (Mi(IV) - l ) 4 + (M 2 (IV) - l ) 4 + 4 
< e (III.154) 


for all N > N, N E Af. Consequently, f nm(n)(jIn) —^ /('?), as IV —* oo, which 
completes the proof of condition (a) of Proposition III. 14. 

In order to show that condition (b) of Proposition III. 14 is satisfied, 
suppose that a sequence {r/^NeAT is such that t/n G for all N E A f, and tjn — r/, 
as At — * tx). Based on the construction of H^, we must have that rj E H” . It again 
follows from the H-continuity of /(•) and Proposition III.26(i) that f nm(n){j]n) —^ 
f(rj) as N —> oo, which satisfies condition (b) of Proposition III.14. This proves that 


(GTP NM ( at)) epi-converges to ( GTP d ). 

For the second part of the proof, we show the convergence of the opti¬ 
mality functions. Suppose that an infinite sequence {t]n} N€ j^ is such that r] N G H° N , 
for all N E A f, and tjn —^ V as N —> oo. Let e > 0. From Lemma III .6 and Lemma 
III.19 we know that V/(-) is Lipschitz H-continuous on bounded subsets of H^. By 
Assumption III. 28 and the H-continuity of V/(-), there exists an N G Af such that 


for a " iiuAH? 


< < j and I|V/(i?k)-V/(j))|| Fs < 


~ — - - - ) (Mi(AT)-l ) 4 1 (M 2 (AT)-1 ) 4 - 2f N - 4 — M v ■' V'uvy * J VU\\n 2 _ 4 , 

where Kn and Kn are as in Proposition III.26, and Kp is as in Proposition III.18. 


Hence, by Proposition 111.26(h), 


\\¥f nm{n){jin) - V/(^)||ff 2 

< II A7f N M(N)(VN) - V/M \\h 2 + ||V/(77jv) - V/(77)|U 2 

K f K n K n e 

“ N + (Mi(IV) - l) 4 + (M 2 (IV) - l) 4 + 4 

< e, (III.155) 


for all N > N, N G A/”. Consequently, V/jVM(JV)(?hv) —V/(r 7 ) as A” —>■ 00 , and 
therefore 0nm(n)(vn) —^ 0{rj) as A” —>• 00 . The epi-convergence of ( GTPnm(n)) to 
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(■ GTP d ) as N —» oo, together with the convergence of Onm(n)(t]n) to 6{rj) as N —» oo, 
satisfy the requirements of Definition III.2 for consistency of approximation, which 
completes the proof. □ 

We now show that the pairs {{GTP^ M , N C),6 c NM , N j) in the sequence 

{((GTF% m{n) ), ^nm(n))}noM are consistent approximations for the pair ((GTP C ), 9 C ). 

Theorem III.30. Suppose that Assumptions III.2, III.3, III. 15, and III.28 are sat¬ 
isfied, ( GTP C ), 6 C , ( GTPf M ^), and 9 C NM ^ are defined as in (III.20), (III.52), 
(III. 151), and (III. 153), respectively. Then the pairs ((GTP^ M , N ^ 9 c NM( ^ N f), in 
the sequence {((GTP^ M , N ^), On M ( N ))}n<=N are consistent approximations for the pair 
(( GTP C ),0 C ). 

Proof. The proof that the problems (GTPfj M , N j) epi-converge to ( GTP C ) is the 
same as the proof that the problems (GTP NM ( N f) epi-converge to ( GTP d ) given in 
Theorem III.29 above. 

From the proof of Theorem III.29 above, we know that V /nm(n) (vn) — 
Vf(r]), as N —> oo. Then, following the same arguments as in the proof of Theorem 
4.3.6 in Polak (1997) we see that given any infinite sequence {vn}n & j^, such that 
rjN € H C;JV for all N G A/", which converges to an r/ G H c , 9 ( )j M ( N )(vn) —^ O c (rf), as 
N —y oo. □ 

C. INDEPENDENT TARGETS 

We now consider the case of independent targets, where we proceed in a man¬ 
ner similar to that used in Section III.B for the case of dependent targets. While our 
assumptions and approach are similar to those used in Section III.B, the “informa¬ 
tion state” for the case of independent targets is different than it was for the case of 
dependent targets. The difference in the “information state” is explained in detail in 
Section III.C.l. 

This section provides discretization schemes that lead to implementable algo¬ 
rithms for the problems ( ITP e ), (JTP c,e ), ( ITP P ), and (JTP C,P ), which were defined 
in Chapter II. The definitions of these problems in Chapter II were incomplete be- 
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cause they did not include definitions for the spaces of allowable controls. This sec¬ 
tion begins by defining an “information state” which we use in conjunction with the 
spaces from Section 111.A to complete the definitions of ( ITP e ), ( ITP c,e ), ( ITP P ), 
and ( ITP C ’ P ). Next, we state our assumptions and define optimality conditions for the 
generalized optimal control problems. Then, we develop consistent approximations 
for the time-discretized search problems. Finally, we show that the time- and space- 
discretized search problems are consistent approximations for the original, continuous 
time-and-spaee search problems. 

1. Information State and Optimal Control Problems 

Under the assumption that the random variables that the target motion is 
conditioned upon are independent across targets, we use (II.9) to define the function 
f l : H —» M for any rj G H by 


K 


f l (v ) = 


exp 


' a l &A 


— / r k ’\x 11 ’ k (t),y l (t', a l ))dt <fi\ot l )da l . (III.156) 


k=\ 


As before, the superscript l is used to denote the I th attacker and l = 1, 2,..., L. Again, 
to simplify the notation in (III. 156) and facilitate the development that follows, we 
find it useful to define a parametric “information state” denoted by z v,l (t ; a 1 ). For any 
a 1 G A, l = 1,2, ...,L, t <E [0,1], and set of searcher trajectories, z 1,,l (t;a l ) represents 
the cumulative detection rate given those searcher trajectories and a 1 , and is given 
by 

r t K 

z v,l (t] a 1 ) = / Vr y (x v ’ k (s),y l (s-,a l ))ds, (III.157) 

Jo k =i 

or equivalently by the differential equation 

I< 

z v,l (s] a 1 ) = J2 rk ’ 1 {x v ’ k {s),y l (s]a 1 )) V s G [0,t], (III.158) 

k=\ 


with z n,l { 0; cd) = 0. It is important to note that the “information state” in (III.157) 
differs from the “information state” in (III.12) for the dependent target case. In the 
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dependent target case, the “information state” represented the cumulative detection 
rate for all of the searchers looking for all of the targets. For the independent target 
case, the “information state” represents the cumulative detection rate for all of the 
searchers looking for the I th target. Using this notation, for any rj G H. (III.156) 
simplifies to 

f l {rj) = f exp (— z v,l (l-, a 1 )) (j) l (a l )da l . (III.159) 

Ja l eA 

Again, it is useful to simplify the notation in (III. 159) even further. To this 
end, for any a 1 G A, l = 1, 2,..., L , we also define the function a 1 ) : H — » M by 


f l (ri;a?)±F(£,xrt(t;o?) J, (III.160) 

where F(-; •) is defined as in (III.15), and x v,l (t ; a 1 ) is an augmented state defined by 

x v (t) 

.z v,l (t ; a 1 ) 

Using this notation, for any r) G H. (III.159) simplifies to 

(III.162) 


x v \t;a l ) = 


G R nI<+1 . 


(III.161) 


f l (v)= I 

J a 1 £ A 

To formulate the problem of maximizing the expected number of targets de¬ 
tected during [0,1], we define the objective function : H —> M for any 77 G H 
by 

L 

>/>'(») =£ f‘M- (HI-163) 

Z =1 

Then, we consider the problem 


(. ITP e ) min ' 0 e ('?/). 

i;£H° 


We also consider the problem 


(III.164) 


(. ITP c ’ e ) min ^ e (r/), (III. 165) 

t)6H c 

which allows for constraints on the control input. It should be noted that the prob¬ 
lems ( ITP e ) and ( ITP c,e ) are both minimization problems despite the fact that they 
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correspond to maximizing the expected number of targets detected. This is because 
the probability that at least one of the searchers detects the I th target during [0,1] is 
given by 1 — f l (r)), and hence the expected number of targets detected during [0,1] 
is given by J2i=i t 1 “ f l (v)] ■ 

To formulate the problem of minimizing the probability that all of the serachers 
fail to detect any of the targets during [0,1], we define the objective function : 
H —> R. for any rj G H by 

L 

r(v) = l[f l (v)- (III.166) 

i=i 

Then, we consider the problems 


(. ITP P ) min ijj p (r]), (III. 167) 

?;SH° 

and 

(. ITP C ’ P ) min ^ p (r/). (III.168) 

rieHc 


2. Optimality Conditions 

In this section, we state our assumptions and give optimality conditions for 
(JTP e ), ( ITP c,e ), ( ITP P ), and ( ITP C,P ). We begin by deriving parameterized differ¬ 
ential equations of the augmented dynamics in terms of the augmented state, x l (t ; a 1 ), 
defined in (III.161). For all / = 1,2, ...,L and t G [0,1] we define 


h l (x(t),u(t ); a 1 ) = 


( h l {x 1 (t ), u 1 (t)) ^ 


h K ( x K (t),u K (t)) 
VEf=i r k \x k {t\y l {t-a l ))J 


e M nA ' +1 . 


(III.169) 


For a given a 1 G A, l = 1,2, ...,L, we write the following parameterized differential 
equation to describe the augmented dynamics 


x\t) = h l (x(t), u(t), a 1 ), t G [0,1], 5^(0) = (III.170) 
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We let x v,l (-) denote the solution of (III.170) when the input is r/ = (£,u), and 
f = (f T , Of. We next state a series of assumptions similar to those given in Section 
III.B.2, beginning with those related to f (•) and y l (-\ •). 

Assumption III.31. We assume that $(■), l = 1, 2 ,L is four times continuously 
differentiable. □ 

Assumption III.32. We assume that Assumption III.2 holds, with a replaced by a 1 , 
for l — 1, 2, L. □ 


The next assumption is related to r k,l (•,.*;) and /r(-,-;-), where we will adopt the 
notation 

/ MfW + 'l n,( + \\T \ 


h x (x(t), u(t); a ) = 


hl(x(t),u(t)) T 

h?(x(t),u(t)) T 


\Efc=i ( xk (t ), y l (f; a 1 ) ) T ) 

where h l x (x(t),u(t); a 1 ) is a (■nK + 1) x n matrix and 


/ 


h l u (x(t), u (t)-,a 


i\ A 


h l u (x{t),u{t)] 

h«(x(t),u(t)) T 

0 


\ 


(III.171) 


(III.172) 


where h l u (x(t),u(t)', a 1 ) is a (nK + 1) x m matrix. 

Assumption III.33. We assume that Assumption III.3 holds with h(-, •; a), h x (•; a), 
and h u (-, •; a) replaced by h l (-, •; a 1 ), h l x (-, •; a 1 ), and h l u (•; a 1 ), respectively. □ 


We next show that if e (-) and f ) p (•) are Gateaux differentiable on H°. 

Proposition III.34. Suppose that Assumptions III.31, III.32, and III.33 are satis¬ 
fied. Then, for any q e H° and Sr] e 

(a) if e (-) has a Gateaux differential Df> e fi 7 ; Sr]) at 77 given by 


h 2 ’ 


Df> e (r]'i Sr]) = (VfAfq),Sq) 


(III.173) 



where the gradient Vfi e (r]) is given by 


Vrmt) = E ( I V a l )(t)<f>\a l )da l ) , Vt G [0,1], (III.174) 

^ \Ja l eA J 


and 


(b) ip p (-) has a Gateaux differential D^ p {j]\ Sr)) at rj given by 

Dfi p (v, &v) = $v)h 2 , 

where the gradient Vfi p (r)) is given by 


(III.175) 


E 

i=i 


l eA 


V r J l (j)\ a l )(t)4>\a l )da l 


n 

j=Mj¥=i 


cPeA 


f^{ry, oP)<j? (o?)da? 


Vt e [ 0 , 1 ]. 


(III.176) 


Proof. Part (a) follows by the same arguments as those used in the proof of Propo¬ 
sition III.5, with /(•;«) replaced by f l (-',a l ), and the recognition that the Gateaux 
derivative of a sum is the sum of the Gateaux derivatives. Part (b) follows in a similar 
manner, with an application of a product rule for Gateaux derivatives. □ 

As in Section III.B.2, we again observe that Proposition III.34 also holds under 
weaker assumptions on <f> l {-), but because we need Assumption III.31 later we adopt 
it here as well. This issue regarding Assumption III.31 also applies elsewhere in this 
chapter. Our next task is to show that Vip e (-) and \70 p (-) are Lipschitz H-continuous 
on bounded subsets of H n . 

Lemma III.35. Suppose that Assumptions III.31, III.32, and III.33 are satisfied, 
then 

(a) the gradient Vi/i e (-) is Lipschitz H-continuous on bounded subsets of H°, and 

(b) the gradient V^ p (-) is Lipschitz H-continuous on bounded subsets of H°. 

Proof. By Lemma III.9, sums and products of Lipschitz H-continuous functions are 
Lipschitz H-continuous on bounded subsets of H". Based on Lemma III.6 we know 
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that the components of (III.174) and (III.176) are all Lipschitz H-continuous, and 
therefore the entire summations in (III.174) and (III.176) are Lipschitz H-continuous. 

□ 

We define nonpositive optimality functions 9 e : H° —>■ M, 0 c,e : H c —> M, 
9 P : H° —>■ M, and 9 C,P : H ( —> M by 

#'(<?) = ~)l|vr(>))llk. (HI-1 77 ) 

8 e ' e W - mill (Vt/i'O), i i - r)) Hl + ) ||i)' - „f„,, (III.178) 

T] tn c Zj 

« p (i?) = (hi. i79) 

and 

9 C ’ P ( V ) = min (V#^), rj - rj) H2 + J || r( - r )||^ a , (III.180) 

Tj tn c Z 

which dehne optimality conditions for ( ITP e ), ( ITP c,e ), ( ITP P ), and ( ITP C,P ), re¬ 
spectively. 

Proposition III. 36. Suppose that Assumptions III.31, III.32, and III.33 are satis¬ 
fied. 

(a) 9 e (-), 9 c,e (-), 9 P (-), and 9 C,P (-) are H°-continuous functions. 

(b) If fj E H° is a local minimizer of ( ITP e ), then 9 e (rj) = 0. 

(c) If rj E H, is a local minimizer of ( ITP c,e ), then d c,e (r)) = 0. 

(d) If rj E H° is a local minimizer of ( ITP P ), then 9 p (fj) = 0. 

(e) If f/ E H c is a local minimizer of ( ITP C,P ), then 9 c,p ('rj ) = 0. 

Proof. The proof follows the same arguments as those for the proof of Theorem 
4.2.3 in Polak (1997), with Proposition III.34 taking the place of Corollary 5.6.9 from 
Polak (1997) and Lemma III.35 taking the place of Theorem 4.1.3 from Polak (1997). 

□ 
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3. Consistent Approximations 

In this section we define the approximating problems ( IT Py), (JTP^ e ), (IT Pfj), 
and (ITP^ p ), and present consistency conditions for them. As in the case of de¬ 
pendent targets, we divide our development into two subsections. Both subsec¬ 
tions develop consistent approximations for the pairs ((/TP e ),d e ), ((/TP c,e ), 6 c,e ), 
(( ITP P ),9 P ), and ((ITP C ’ P ),9 C ' P ), but the first subsection only deals with time dis¬ 
cretization while the second subsection considers time and space discretization. 

a. Time-Discretized Problems 

We again consider the approximate solution of (11.22) by means of for¬ 
ward Euler’s method, which was given in (III.60). Simultaneously, we approximately 
solve (III. 158) also by forward Euler’s method. For any r] = (£,w) G Hjy, a 1 G A, and 
N G A f, we set 0; a 1 ) — 0, l — 1, 2,..., L, and for any j = 0,1,..., N — 1, 

1 K 

(U + 1)/A (j/N;a 1 ) = = £>*’' ^(j/N),y l (j/N-,a 1 )) . (III. 181) 

V k =1 


Using the discretized “information state” given by the recursion (III. 181), 


we define the approximate objective functions ^ e,N : 

for any rj G Hat and N £ AT by 

Hat —y M and ij} p,N 

\ —y M 

^n(v) ~ l exp ( _ ^v ( 1; a *< 
1=1 V 

( f> l {a l )da l , 

(III.182) 

and 

L 



v N (v )— n / exp ai ) 

i=i j ^ a v 

4> l (a l )da l . 

(III.183) 

We also define the individual components of (III.182) and (III.183) by 



f N (v) = J ^exp <f> l (a l )da l yi = 1 , 2 , (III.184) 

Again, for the sake of notational simplification, for any a 1 G A, l — 1, 2,..., L, we also 
define the functions f l N (-; a 1 ) : D° N — > M by 

P N (v,a l )^F{i,x^(f,a 1 )), (III.185) 
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where F is as defined in (III. 15) and x v J'(j/N] a 1 ) is an augmented state defined by 


^(j/N-a 1 ) = 


xUj/n) 


L , j — 0,1 ,N — 1. 


W / ' / I 7 , I I III 

YnU/N- a 1 )) 

Hence, for any N G A/", we define the following approximating problems 


(III.186) 



min 

04(h), 

(III.187) 

(upy) 

min 

V^c,N 

04(h), 

(III.188) 

(UPJ) 

min 

04(h), 

(III.189) 

(/TP^) 

min 

H c .jv 

04 (h)- 

(III.190) 


We next consider the differentiability of if ^(•) and 

Proposition III.37. Suppose that Assumptions III.31, III.32, and III.33 are satis¬ 
fied, and N G A I. Then, for any rj G and dp G h/ 00,2 

(a) if%(-) has a Gateaux differential DiJf N (p] dp) = {Vif%(p), dp) H2 , where 

Vr N (v)(t) = E [ ^A^a l )(t)fi\a l )da l yt G [0,1], (III.191) 

1=1 


(b) if p N {-) has a Gateaux differential Dif p N (p\ dp) = {Shf v N fp), dp) Hn , where 

Vl&fa)(f) 

= ( / V r,f l N(V',<x l )(t)(l> l (ai l )da l JJ f f N (p;a J )fi , (a J )da J 

1=1 \ Ua l eA J J = 1 | 


Vt G [0,1]. 


(III.192) 
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Proof. Part (a) can be proven by the same arguments used in the proof of Propo¬ 
sition III. 5 , with /(•;«) replaced by f l N (-]a l ), and the recognition that the Gateaux 
derivative of a sum is the sum of the Gateaux derivatives. Part (b) follows in a similar 
manner, with an application of a product rule for Gateaux derivatives. □ 

Next we show that and V^(-) are Lipschitz H^-continuous 

on bounded subsets of H 1 ^. 

Lemma III.38. Suppose that Assumptions III.31, III.32, and III.33 are satisfied, 
then 

(a) the gradient Wfi> e N {-) is Lipschitz Hjv-continuous on bounded subsets of H° N , and 

(b) the gradient V fi p N (•) is Lipschitz H.y-continuous on bounded subsets of H^. 

Proof. By Lemma III.9, sums and products of Lipschitz Hy-contirmous functions are 
Lipschitz Hy-continuous on bounded subsets of H”,. Based on Lemma III.10(b) we 
know that the components of (III.191) and (III.192) are all Lipschitz Hjv-continuous, 
and therefore the entire summations in (III.191) and (III.192) are Lipschitz H,y- 


continuous. □ 

For any N e Af, we define nonpositive optimality functions 6% : H ( y —* 
R, 6% e : H cJV -)■ R, 0 P N : U° N —)■ R, and Ofifi : H ciV -> R by 

0%(v) = -i||V*0)llk. ( IIL193 ) 

««'(>))- min + (III.194) 

e r M = -i||v*(>))llk. (ni-iss) 

and 

° C n (v) - mj.n W P M,V' - v) h 2 + \ \W - v\\ h 2 , (III.196) 

T) Grl c ,JV £ 


which characterize stationary points of (IT Pfi), (IT Pfifi), ( ITP ^), and {ITPfifi), 
respectively. 
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Proposition III. 39. Suppose that Assumptions III. 31, III.32, arid III.33 are satis¬ 
fied. 

(a) 9 e N (•), 9% e (-), 9 P N (-), and 9 c fif(-) are H^-continuous functions. 

(b) If fj G is a local minimizer of ( ITP ^), then °n(v) = 0. 

(c) If f] G H Ci jv is a local minimizer of (ITPffi), then 9ff(rj) = 0. 

(d) If fj G U° N is a local minimizer of ( ITP’jf), then 9 p N (fi) = 0. 

(e) If f) G H Cj jv is a local minimizer of (ITPff), then 9 c ffi(rj) = 0. 

Proof. The proof follows the same arguments as the proof of Proposition 1.1.6 in 
Polak (1997), with the norms and inner products replaced with their H 2 equivalents. 

□ 

As is the case in Section III.B.3a with epi-convergence of (GTP^) to 
(<GTP ), it is not possible to establish epi-convergence of ( ITP%) to ( ITP S ) or (IT Pfi) 
to (ITP P ). We define the problems 


(ITPS) 

min 

^(V), 

(III.197) 

r?£H0 ; 

(ITPS 

min 

V(v), 

(III.198) 

-)£H ° cl 


because it is possible to establish epi-convergence of ( ITP%) to (. ITP and ( ITP'^) 
to (ITPff). We will use the problems (IT Pf) and (IT Pfi), with the following assump¬ 
tion. 

Assumption III.40. We assume that all local and global minimizers of (ITPfi) and 
(ITP p cl ) are m H°. □ 

Again, in a manner similar to that of Section 3.3 of Polak (1997) we 
next show that the pairs ((ITPff) , 0 e N ) in the sequence {((ITPffi), 9%)}n<=n and the 
pairs ((ITPfi), 9 P N ) in the sequence {((ITPff), 9%)}NeM are consistent approximations 
for the pairs ((ITPff), 9 e ) and ((ITPff), 9 P ), respectively. In order to establish epi- 
convergence, we will need the following intermediate result. 
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Lemma III. 41. Suppose that f3 l : H —y [0,1] and B l N : —* [0,1], l = 1,2, 

N eN, are such that for all rj E H, l — 1, 2,..., L, and N G N, 

\P 1 N (V)~P 1 (V)\<* (N) (III.199) 

where A : N —> [0, oo) is a strictly decreasing function with A (N) —* 0, as N —* oo. 
T/z-en for all N E N and 77 G H, 

(a) 

L L 

/3 Z ( 77 ) < LA(N) (III.200) 

1=1 1=1 

and 

(b) there is a constant K n < 00 such that 

L L 

nAori-n P l (r])\<K v A{N). (III.201) 

1=1 1=1 

Proof. To prove part (a), we expand the summations in (III.200), combine terms, 
and use the triangle inequality to write 

L L L 

= E (Aw - Pw) 

i=i 1=1 1=1 

L 

< 

1=1 

< LA(N) (III.202) 

For the proof of part (b), we write out both the products in (III.201) 
and express /3 l N (r]) in terms of j3 l (r]) and A(IV), then because /3\r]) E [0,1] and 
/ 3 l N (r) ) G [ 0 , 1 ] for all l we have that 

L L 

nAM-rpM 
1=1 1=1 



(111.203) 

(111.204) 
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where the final inequality follows from taking the product of the (/ 3 l (rj ) + A (TV)) in 
(III.203) and combining like terms. Then because A (TV) —>■ 0 as N —>■ oo, there exists 
N e N such that for all N > N there is a constant K n < oo such that 

E (f) = “W + E (t) ^ A'.A(iV), (III. 205 ) 

Z=1 ' ' 1=2 ' ' 

□ 

We now show the epi-convergence of ( ITP to ( ITP and ( ITP 

to (ITPS). 

Theorem III.42. Suppose that Assumptions III.31, III.32, III.33, and III.40 are 
satisfied. Then 

(a) (ITPfj) epi-converges to (ITPff), as N —> oo, and 

(b) {ITPffi) epi-converges to (ITPfi), as N —> oo. 

Proof. The proof follows the same arguments as the proof of Theorem III.17. We 
then invoke Lemma 111.41(a) and (b) to complete the proofs of parts (a) and (b), 
respectively. 


□ 

In order to show consistency of approximation for the pairs ((/TP^), 9 e N ), 
(filTPff), 9ff), (( ITPfj ), 9 p n ), and ((ITPff), Off) we will need the following two 
results. 

Lemma III. 43. Suppose that there exist constants, C l < oo and C l N < oo, and that 
the functions V: H° —> H. V f3 l N : H^r —> H^v, : H —> [0,1], and (3 l N : —> 

[0,1], l — 1, 2,...... L, N e N, satisfy 

||V /^)|\ H2 < C l ,Vl = l,2,...,L, V ett° (HI-206) 

\\VPM\\h 2 < C‘ ) VPl ) 2 1 .,L, f) 6H»,iVeN (III.207) 

\Pn(v)~P 1 (v)\ < A(IV), \/l — 1, 2,..., L,rj £ Hjv, IV G N (III.208) 

and 

\\V/3 l N (v) - V/3 l (r])\\ H2 < A(IV), V/ = 1,2,..., L, 77 e H ° N ,Ne N, (IIL209) 

where A : Id —> [0, 00 ) is a strictly decreasing function with A(IV) —> 0, as N —> 00 . 
Then, 
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(a) for all N G N and r) G 


L 


L 


E v &M-E v rt?) 


< LV2A(N), 

h 2 


and 

(b) there is a constant K^g < oo such that for all N G Id and rj G 


(III.210) 


e v aw n n 


i=i 


j=i|jyy 


z=i 




< K vG A(N). 


h 2 


(III.211) 


Proof. We begin by proving part (a). We know that 




1=1 


i=i 


L 

<EII v A(>))-v/j , o)|| H3 . 

H 2 1=1 


(III.212) 


Let / = 1,2,..., L, IV G N, and r) G be arbitrary, then based on the definition of 


the H 2 norm || V/3 l N (r]) — ^P l {v)\\ H2 < A(IV) implies that 

||vA(r,)-V/3'(,)|| 2 Hi = ||V £ A(ij) - V 5( S'(i))|| 2 

+ / IIWAWW-V^WWII 2 * 

do 

< [A(IV )] 2 + [A(IV )] 2 , (III.213) 


where S l N (r /) G M n , V^f3 l (r]) G M n , V u (3 l N (r]) G U, and V u /3 l (r)) G U. Hence 

||V^(77) - V/3 i (7 7 )|| jff2 < V^A(IV),VI = 1,2,..., L, (III.214) 

which completes the proof of part (a). 

For the proof of part (b), let g l N {rj) A /%(?l) and ^('7) - 

YlUwFW- Then ’ 

L L L L 

evaw n n 

i = 1 l=lbW i = 1 rr 
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= E - Y s'wvp'w 

1=1 1=1 h 2 

L 

< Yb'MVpM-g'WVP'MWK- (HI-215) 

(=1 

Based on the definition of the H -2 norm we have 

IIAMVAO) -s'OW'Mllh = 

||9t(>f)V £ A0) - s'WV^'WII 2 + f IIAWVAWW - g'(vW(ri)(t)\\ 2 dt 

Jo 

(III.216) 

Based on Lemma 111.41(b), there exist constants K l < oo, l = 1,2 , ...,L, such that 
V7V e N, and Vr? e 

WnM - g‘(ri)\ < K‘A(N). (III.217) 

We then consider the first term on the right-hand side of (III.216), which can be 
rewritten 

= IkMvAM - g'MVtP'W + (g'M - g‘M) v £ /?U>i)|| 2 

< g'lnf ||v £ Ah) - v £ /3‘(,)|| 2 + ||(<4(>?) - g'w) v £ A (,)|| 2 

< [A(AQ] 2 + (gUv) - g‘M) 2 1| v £ AW|| 2 

< [A(AI)1 2 + (Kg) 2 |A(/V)] 2 IIVjAWII 2 

< [A(AI)] 2 + (Cjy) 2 (Kg) 2 [A(A7)] 2 

< (l + (C-i-) 2 (A'J) 2 ) [A(7V)] 2 . (III.218) 
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Next, we consider the second term in the sum on the right-hand side 
of (III.216), which can be rewritten 

/ IKMvAMW - g'inWuP'^mfdt 

J 0 

= [ || g l (v)v u P l N (v)(t) - g l (v)v u P l (v)(t) + ( g l N (v) - g\v)) v u p l N (v)(t)\\ 2 dt 

Jo 

< g‘(v) 2 [ ||v„AW(t) —v„/3‘(r,)(i)|| 2 di+ / IKAW-s'W) vAinmfdt 

Jo Jo 

< [A(JV)] 2 + f\g‘M-g‘M) 2 \\vA(vm\\ 2 dt 

Jo 

< [A(/V)] 2 +(A'J) 2 [A(Af)] 2 f'\\V u P‘M(t)\\ 2 dt 

Jo 

< [A(JV)] 2 +((C',) 2 (Aj) 2 )[A(Af)] 2 < (l + (C‘,) 2 (Aj) 2 )[A(Af)] 2 . (III.219) 

Hence 

IIAMVAM - g'WVtfWW^ < ^2 (l + (Cf ,) 2 (Ag) 2 ) A(iV), V/ = 1,2, L. 

(III.220) 

Then K n c = 22t=i \Jz ^1 + (Cat)" (A^J, and the proof is complete. □ 

Proposition III.44. Suppose that Assumptions III.31, III.32, and III.33 are satis¬ 
fied, then for every bounded subset S C H°, 

(a) there exists a constant K^p < oo such that, for any N e M and rj G S ft H N 

\\vrM-vrm\H, < 'ffi ( 111 . 221 ) 

and 

(b) there exists a constant K n p < 00 such that, for any N e J\f and 77 £ S fl Hn 

II - VV p (»)|k < ^. ( 111 . 222 ) 

Proof. From Proposition III. 16, we deduce that 

[ Pn(v, a 1 ) - f l (m a 1 ) </>VW < ^ • (IH.223) 

da'eu iV 
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The proof then follows the same arguments as the proof of Proposition III. 18 with 
Propositions III.5 and III.11 replaced by Propositions III.34 and III.37, respectively. 
We then invoke Lemma 111.43(a) and (b) to complete the proofs of parts (a) and (b), 
respectively. □ 

We now show that the pairs (( ITP%) ) @n) in the sequence 
{{{ITPff),9 e N )} Ne u and the pairs {{ITP^),9 p N ) in the sequence {{{ITP^),9 p N )} Ne jy 
are consistent approximations for the pairs {{ITP^),9 e ) and (( ITP^),9 P ), respec¬ 
tively. 

Theorem III. 45. Suppose that Assumptions III.31, III.32, III.33, and III.40 are 
satisfied, ( ITP%), ( ITP p ), 9 e , 9 P , {ITPffi), ( ITP'^), 9 e N , and 9 P N are defined as in 
(III. 197), (III. 198), (III. 177), (III. 179), (III. 187), (III. 189), (III 193), and (III. 195), 
respectively. Then 

(a) the pairs {{ITPffi), 9 e N ), in the sequence {((/TP^),^)}^ are consistent ap¬ 
proximations for the pair {{ITP^fi, 9 e ), and 

(b) the pairs {{ITP((), 9 P N ), in the sequence {((ITPjy), 9 p n )}n&M are consistent ap¬ 
proximations for the pair {{ITP p cl), 9 P ). 

Proof. The proof follows the same arguments as the proof of Theorem III.20 with 
Theorem III.17, Lemma III. 6 , and Proposition III.18 replaced by Theorem III.42, 
Lemma III.35, and Proposition III.44, respectively. □ 

We next show that the pairs {{IT Pfi e ) ^) in the sequence 
{{{ITP c N e ),9 c ^)} N£Af and the pairs {{ITPfifi), 9 c fi p ) in the sequence {{{ITP^- p ), 9 c fifi)} NgJ ^ 
are consistent approximations for the pairs (( ITP c,e ),9 c,e ) and {{ITP C,P ), 9 C,P ), respec¬ 
tively. 

Theorem III. 46. Suppose that Assumptions III.31, III.32, III.33, and III.40 are sat¬ 
isfied, {.ITP c ’ e ), {ITP C,P ), 9 c ' e , 9 C ’ P , ( ITP)) e ), {ITPfifi), 9% e , and9% p are defined as m 
(III. 165), (III. 168), (III. 178), (III. 180), (III. 188), (III. 190), (III. 194), arid (HI-196), 
respectively. Then 

(a) the pairs {{ITP^ e ), 6^ e ), in the sequence {{{ITP^ e ), are consistent 

approximations for the pair (( ITP c,e ),9 c ' e ), and 

(b) the pairs {{ITP^f), 9 c jf), in the sequence {{{ITPfif),9 c I (fi)} are consistent 
approximations for the pair {{ITP C,P ),9 C,P ). 
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Proof. The proof follows the same arguments as the proof of Theorem III.21 with 
Theorems III. 17 and III.20 replaced by Theorems III.42 and III.45, respectively. □ 


b. Time- and Space-Discretized Problems 

We next consider the time- and space-discretized problem. We make 
use of the integration rule Im given in (III. 102) to define the approximate objective 
functions : H^v —» M and fi P NM : —> M for any rj G Hn, N e Af , and 

Me NxNby 


^nm(v) - Im ( exp ~ z n\ 1 \ •) , 


and 


i=i 

L 


vnm(v) - n Im ( exp - z n i ( v ’ 


1=1 


(III.224) 


(III.225) 


We next consider the differentiability of and V ; atm(')- 

Proposition III.47. Suppose that Assumptions III.31, III.32, and III.33 are satis¬ 
fied, N e Af, and M e N x N, then for any r) e and Sr) e 

(a) h as a Gateaux differential Dfi e NM (r)\ Sr)) = ^)# 2 , w h ere 


= £ (Im [v,A(d;-)W^(-) 1 ) ,vt e [o.i], (m.226) 


1=1 


and 


(b) V’atm(') Pas a Gateaux differential Sr)) = (V^ M (r/ ),Srj) H , where 


Vfe(#) = X) 

z=i 



L 

I M [vJ l N (r)r)(t)<i> l (-)W 

n iM fN(v,w(-) 


_l=i \j¥=i 


Vt e [0,1]. 


(III.227) 


Proof. The proof of part (a) follows the same arguments used in the proof of Propo¬ 
sition III.22, with /(•;«) replaced by f\ y(-;cd), and with the recognition that the 
Gateaux derivative of a sum is the sum of the Gateaux derivatives. Part (b) follows 
in a similar manner, with an application of a product rule for Gateaux derivatives. □ 
Next, we show that V-z/’a tm(') an( l ^ 7 '0atm(') are Lipschitz H^-continuous 
on bounded subsets of H° N . 
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Lemma III.48. Suppose that Assumptions III.31, III.32, and III.33 are satisfied, 
then 

(a) the gradient V^m^) is Lipschitz H^-continuous on bounded subsets of H° N , 
and 

(b) the gradient V^ M (?/) is Lipschitz H^r-continuous on bounded subsets of H° N . 

Proof. The proof follows the same arguments as the proof of Lemma III.38, with 
integration replaced by Im- □ 

For the reasons discussed in Section III.B.3b, in order to conduct the 
analysis that follows we again make the assumptions delineated in Assumption III.24. 
Based on the choice of Composite Simpson’s rule as the numerical integration scheme, 
we find it necessary to show that the partial derivatives of f l N (ry, ■)4 > l (•) and 
VMm ')^(') U P to and including the fourth-order are bounded for any choice of 
r) G any a 1 G A, and any N G A f. The bounding of the partial derivatives 
of f l N (ip, -)<f> l (-) and V T ,/^(r/; -)0 Z (-) is necessary in order to prove the convergence of 
^nm(v) to ^(v), Vfi e NM (v) to Vr(v), ip P NM (V ) to and Vfi p NM (r)) to Vfi p (r]). 

To facilitate these proofs we begin by defining some notation. For any rj G H%, 
a 1 G A, l = 1, 2,..., L, N G Af ', and j = 0,1,..., N — 1, we define 


C1K) ^ exp 





(III.228) 




(III.229) 


Ci( al ) — Pjv (0> cv 1 ), (III.230) 

and 

<>") = Pn <*'). (III.231) 

We again note that (({(•), (l('), and Ci(')i ^ = 1,2, ...,L, depend on rj and N. 
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We next show that the partial derivatives of ([ up to and including the 

fourth-order are continuous and bounded for any choice of rj G a 1 G A, l = 
1, 2, L, and N G A f. 

Lemma III. 49. Suppose that Assumptions III.31, III.32, and III.33 are satisfied and 
S is a bounded subset ofH%. Then, 

(i) C'(-)eC 4 (i) Vi = 1,2, 3,4, 
and 

(ii) there exists C < oo, such that for all r/ G S, j — 0,1, ...,N — 1, a 1 G A, 
l — 1,2,L, and N G A/” 

< C Vi = 1, 2, V/i = 1, 2,3,4, V/c = 1,2, 3,4. (III.232) 

Proof. The proof follows the same arguments as the proof of Lemma III.25 with 
Assumptions III. 1, III.2, and III.3 replaced by Assumptions III.31, III.32, and III.33, 
respectively. □ 

In order to prove epi-convergence and consistency of approximation, we 
will need the following proposition. 

Proposition III. 50. Suppose that Assumptions III.24, III.31, III.32, and III.33 are 
satisfied. Then for every bounded subset S C H°, there exist constants Kes < oo, 
K n s < oo, and K n j\ , K n j 2 < oo such that, for any N G A f, for any M G N 3 x Ns 4 , 
and q G S fl H N , 

(a) (i) | ^ e (v) ~ VnmW )| < ^ r + 

(ii) || V^ e (r/) - VVnMWh 2 < ^ + + pifwV’ 

and 

(b) (i) |V> p (ij) - V>S™MI < V + (jIV 1 + (W 1 ’ 

(ii) ||V^(,,) - V«„(r;)|| Hj < Y + wA? + ffiSv 
and K-ef and K n p are defined as in Proposition III. 44- 

4 Recall that N 3 = {m € 2N + l|m > 3}, as defined previously in Proposition III.26. 


^C(a l ) 

dolf 
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Proof. The proof follows the same arguments as the proof of Proposition III.26 with 
Proposition III.18, Lemma III.25, and (i(-),i^-) replaced by Proposition III.44, 
Lemma III.49, and Ci(')> - - -, Cl (*) j respectively. □ 

For any N e AT, and MeNxN, we define the following approximating 

problems 


(ITP e NM ) 

min 

V£K° n 


(III. 233) 

(ITPTm) 

min 

i)6H CiW 


(III. 234) 

( ITP% m ) 

min 

*?6 hS, 


(III. 235) 

and 




(ITP c / m ) 

min 

»?eH CiJV 

y )• 

(III. 236) 

For any N e AT, and M 6 N x N, we define nonpositive optimality 

functions 9 e NM : -x M, 9% e M : H c N 

— x M, 9 p nm \ H^r —x K., and 9 

JYM : H C ,JV —X M 

by 




@nm( t , i) — ■ 

-\w 

Nm(v) \\h 2 i 

(III. 237) 


0°nm(v) = nun (Vip e NM (r]), rf - rj) H2 + \ II V* - v\\ h 2 > (III.238) 

»««(>)) = -l\W NM m 2 H„ (III.239) 

and 

9 n P m(v) = nun ~ v)h 2 + \ II v' - v\\ h 2 , (III.240) 

PeH c ,n 2 

which characterize stationary points of (ITP^ M ), (ITP^ e M ), (ITP^ M ), and (ITP^ P M ), 
respectively. 
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Proposition III.51. Suppose that Assumptions III.24, III.32, arid III.33 are satis¬ 
fied. 

(a) 9 e NM (•), 9 c fi e M (-), 9 P NM ( •), and 9 c fifi M (-) arc H^-continuous functions. 

(b) If fj G is a local minimizer of (ITP^ M ), then 9 e NM (f)) = 0. 

(c) If f) G H c at is a local minimizer of (ITP^ e M ), then @nm(v) ~ 0- 

(d) If p G is a local minimizer of (ITP^ M ), then 9 p NM (rj) = 0. 

(e) If f) G H c at is a local minimizer of (ITP^- P M ), then = o. 

Proof. The proof follows the same arguments as the proof of Proposition 1.1.6 in 
Polak (1997), with the norms and inner products replaced by their H -2 equivalents. □ 
It is not possible to establish epi-convergence of (ITP^ M ^) to ( ITP e ) 
or (ITP^ ) M (at)) ( ITP P ), but it is possible to establish epi-convergence of (ITP^ M ^) 

to ( ITP and (ITP(( M ^) to ( ITP ( p ). We now show that the family 
{((ITPn M ( N )) >^m()v))K a/" is a sequence of consistent approximations for 
((ITP ^), 6 e ), and the family {((ITP^ m ^ n ^), 9 p nm ^)}ngM is a sequence of consistent 
approximations for (( ITP p ), 9 P ) by showing that the conditions of Definition III.2 are 
satisfied. 

Theorem III.52. Suppose that Assumptions III.28, III.32, III.33, and III.40 are sat¬ 
isfied, (ITP*), (ITPPcl), 9 e , 9 P , ( ITP e NM{N) ), (ITP p nm{n) ), 9 e NM(N) , and 9 P NM(N) are 
defined as in (III. 197), (III. 198), (III. 177), (III. 179), (III.233), (III.235), (III.237), 
and (III. 239), respectively. Then 

(a) the pairs ((ITP^ M{N] ), 9 e NM{N) ), in the sequence {((ITP^ M{N) ),9 e NM{N) )} Ne ^ 
are consistent approximations for the pair ((ITP^),9 e ), and 

(b) the pairs ((ITP p NM{N) ), 9 P NM{N) ), in the sequence {((/TP£ M(JV) ), 9 p NM{N) )} NeAf 
are consistent approximations for the pair ((ITP p (), 9 P ). 

Proof. The proof follows the same arguments as the proof of Theorem III.29 with 
Proposition III.26 and Theorems III. 17 and III.20 replaced by Proposition III.50 and 
Theorems III.42 and III.45, respectively. □ 

We now show that the pairs ((ITP^ e M , N ^),9 c ^ M , N ^) in the sequence 
WTP c N ‘ Mm ) , ^v'M(V))} Ar e M are consistent approximations for the pair ((PTP c,e ), 9 c,e ). 
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We also show that the pairs {{ITP^ p m , n ^9 C j^ m ^ n ^) in the sequence 
i((ITP c / M{N) ), 0 C nm(n))} n £M are consistent approximations for the pair ((PTP C,P ), 9 C,P ). 


III. 33, and III .40 are sat- 
9 


Theorem III.53. Suppose that Assumptions III. 28, III.I 
isfied, ( ITP «), (ITP C - P ), 0“ <?*, (ITP" Mm ), ( ITP‘ NM(N)I 

’ U NM(N)> 

are defined as in (III.165), (III.168), (III.178), (III.180), (III.234), (HI-236), (III.238), 
and (III.240), respectively. Then 


and 9 C,P 
ana v NM ( N ) 


(a) the pairs ((/TP^ (JV) ), 9 c N e M{N) ), in the sequence {((ITP^ e M{N) ), 9 c £ m{n) )} n& m 
are consistent approximations for the pair ((/TP c,e ), 9 c,e ), and 


(b) the pairs (( ITP C / M{N) ), 0 C * M{N) ), in the sequence {({ITP c / M{N) ),9 c ^ M{N) )} NeH 
are consistent approximations for the pair (( ITP C,P ),9 C,P ). 


Proof. The proof follows the same arguments as the proof of Theorem III.30, with 
Theorem III.29 replaced by Theorem III.52. □ 



IV. 


RATE OF CONVERGENCE ANALYSIS 


A. INTRODUCTION 

In this chapter, we develop rate of convergence results by expressing the rate 
of convergence in terms of the computational work, rather than the typical number of 
iterations or level of discretization. We relate the compuational work to the number 
of iterations as well as to the levels of discretization by making computational work 
assumptions. This relation allows us to determine a guaranteed rate of convergence 
for the error between the objective function evaluated at iterates generated from the 
discretized problems and the optimal value of the original problem as a function of the 
computational work. Because these results are applicable to a range of approximation 
problems, we prefer to develop our results using abstract problems as we did in Section 
III.B.3. We again define the infinite-dimensional problem 

(P) min f(x), (IV. 1) 

x£X 

where /(•) and X are defined as in Section III.B.3a. For all N, M e N, let / nm '■ 
Bn — * M be a lower semicontinuous function that approximates /(•) on Bn , and let 
X N C Bn be an approximation to X, where Bn and M are defined as in Section 
III.B.3a. We then define the family of finite dimensional approximating problems 

(Pnm) min f NM (x), N, M e J\f. (IV.2) 

x£Xn 

We note that x is now a decision variable, and it no longer represents the physical state 
of the searcher. We also note that the problems (GTP), ( GTP C ), (FTP 6 ), (JTP c,e ), 
(ITP P ), and (ITP C,P ) discussed in Chapter III are examples of (P). The problems 
(GTP nm ),(GTP c nm ), ( ITP e NM ), ( ITP c / m ), ( ITP p NM ), and (. ITP C / M ) discussed in 
Chapter III are examples of (Pnm) under the assumption that Mi = M 2 = M. 
Although the examples we provide from Chapter III are all generalized optimal control 
problems, it is worth noting that the results of this chapter are applicable to other 
types of problems as well. The results may be applicable to approximation problems 
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for (P) that involve a single discretization parameter for the space of decision variables 
and two discretization parameters for the objective function. 

To solve the infinite-dimensional problem (P) some form of discretization is 
necessary. If. for example, the infinite-dimensional problem takes the form of ( GTP) 
and N and M represent time and space discretization levels, respectively, then the 
discretization in time of the searcher’s dynamics using Euler’s method and approxima¬ 
tion of spatial integration via Simpson’s rule described in Section III.B.3b combined 
with a standard nonlinear programming algorithm could be used to find a solution 
to (P). We refer to this type of solution methodology as a discretization algorithm. 

In order to solve (P) using a discretization algorithm, there is a fundamental 
trade-off between the level of discretization and the computational work required to 
approximately solve the resulting finite-dimensional problem ( Pjvm )• A fine level of 
N and M discretization ensures that ( P/vm ) approximates (P), in some sense, with 
a high degree of accuracy. The issue with this approach is that the function and 
gradient evaluations in ( Pnm ) tend to be expensive, and therefore the computational 
work required is high. A lower level of N and M discretization results in a faster 
solution time, due to the relatively less expensive function and gradient evaluations, 
but at the expense of a less accurate approximation of (P). It is usually difficult in 
practice to effectively manage the trade-off between these diametrically opposed goals 
of solution accuracy and computational cost. 

In this chapter, we investigate the rate of convergence of a class of discretiza¬ 
tion algorithms as a computing budget tends to infinity. We show that the policy for 
selecting N and M discretization levels relative to the size of the available computa¬ 
tional budget influences the rate of convergence of discretization algorithms. We iden¬ 
tify optimal discretization policies, in a specifically defined sense, for discretization 
algorithms used to solve the resulting finite-dimensional problems based on finitely, 
superlinearly, linearly, and sublinearly convergent optimization algorithms. 

Other researchers have examined the rate of convergence for optimization prob- 



lems where discretization is necessary to obtain solutions. Dunn and Sachs (1983) 
considers the effect of approximate function and gradient evaluations on the con¬ 
vergence rates of the conditional gradient and projected gradient optimization algo¬ 
rithms. Kelley and Sachs (1986) examines the effect of approximating the Hessian on 
the convergence rate of the BFGS-method when it is used to solve discretized opti¬ 
mal control problems. Dupuis and James (1998) presents a method for obtaining rate 
of convergence estimates for finite difference approximation schemes for stochastic 
and deterministic optimal control problems. The method given in Dupuis and James 
(1998) requires some assumption about the smoothness of the objective function be¬ 
ing approximated, which they exploit to obtain sharp rates and in some cases an 
expansion of the discretization error in terms of the discretization parameter. Kang 
(2008) gives rate of convergence results for the approximate optimal cost computed 
when using pseudospectral methods to solve optimal control problems. These studies 
all share two important similarities. The first being that the approximation schemes 
under consideration are all based on a single discretization parameter, and the sec¬ 
ond being that none of the studies consider the error introduced by the optimization 
algorithm in their determination of the rate of convergence. 

More recent studies in the Monte Carlo simulation and simulation optimization 
literature (see Chen & Shi, 2008 for a summary) deal with how to optimally allocate 
a computational budget across different tasks within the simulation and to determine 
the resulting rate of convergence of an estimator as the computational budget tends to 
infinity. In the stochastic optimization literature, Pasupathy (2010) and Royset and 
Szechtman (2011) consider optimization algorithms with sublinear, linear, and super- 
linear rates of convergence for use with the sample average approximation approach, 
determine optimal policies for allocating a compuational budget between sampling 
and optimization, and quantify the associated rate of convergence of the sample av¬ 
erage approximation approach as the computational budget tends to infinity. In the 
semi-infinite minimax literature, Royset and Pee (2011) examines the convergence 



rate as the computing budget tends to infinity, and provide allocation polices that 
maximize the convergence rate of sublinear, linear, superlinear, and finite algorithms 
used to solve semi-infinite minimax problems. Royset and Pee (2011) also examines 
the solution of the semi-infinite minimax problem by exponential smoothing algo¬ 
rithms, provides an optimal discretization and smoothing policy, and determines the 
corresponding rate of convergence as the computing budget tends to infinity. We note 
that the section on smoothing algorithms in Royset and Pee (2011) uses two distinct 
discretization parameters. Our approach in this chapter follows the methodology of 
Royset and Szechtman (2011) and Royset and Pee (2011). Our analysis is based on 
different assumptions, however, and as a result we reach different conclusions. 

The organization of this chapter is as follows. We begin by defining terminol¬ 
ogy and presenting our assumptions. Next we consider finite, superlinear, linear, and 
sublinear optimization algorithms for solving the finite-dimensional problem ( Pnm )• 
For each class of optimization algorithm we give an optimal discretization policy with 
corresponding rate of convergence, as the computational budget tends to infinity. 
Finally, we state our conclusions. 


B. TERMINOLOGY AND ASSUMPTIONS 


We begin with some assumptions about (P) and (Pnm)- 
Assumption IV. 1. We assume that the following hold: 

(i) The set of optimal solutions X* of (P) is nonempty. 

(ii) There exist constants N, M G N, and K G [0, oo) such that 

(a) the set of optimal solutions X^ M of (Pnm) is nonempty for all N > N, 
M > M, A, M G N, and 

(b) there exist p, q G (0, oo) such that 


\f(x) ~ f NM (x)\ < ^ + 


(IV.3) 


for all x G X N , N > N, M > M, and N, M G N. 
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If, for example, ( Pnm ) is a generalized optimal control problem of the form defined in 
Chapter III, then part (b) of item (ii) holds when (Pnm) is solved using an algorithm 
that utilizes Euler’s method to approximate the solutions of the differential equations 
and approximates spatial integration using Simpson’s rule. In this case p = 1 and 
q = 4; see Proposition III.26 and Proposition III.50. We refer to 

f(x) - f nm(x ) (IV.4) 

as the discretization error. Unless Xn and / nm(x ), x G Xn, have special structures, it 
is impossible to obtain a globally optimal solution of (Pnm) in finite computing time. 
Therefore, once a finite number of iterations of an optimization algorithm are applied 
to (Pnm), there is usually optimization error. Given an optimization algorithm A 
for (Pnm), let x^ M 6 Xjv be the iterate obtained by A following n iterations on 
(Pnm)- We use the notation to denote the optimal value of (Pnm) and define 
the optimization error as 

fN m(x^n M ) ~ fNM' ( IV . 5 ) 

The rate of decay of the optimization error as n gets larger depends on the rate of 

convergence of A when used to solve (Pnm)- Throughout this dissertation we only 
consider deterministic algorithms that generate iterates in A" exclusively, which we 
formalize in the next assumption. 

Assumption IV.2. We assume that if {x nm } ti = o > N, M G N, are generated by a 
given optimization algorithm when applied to (Pnm), then x^n G X for all N,M G N 
and n = 0,1, 2 ,.... □ 

Based on the definition of X in Section III.B.3a, many optimization algorithms 
satisfy Assumption IV.2. We use the notation f* to denote the optimal value of (P) 
and define the total error as 

I/Km)-A (IV.6) 
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which is a measure of the quality of the solution obtained after applying n iterations 
of A to ( Pnm )• Based on Assumptions IV.1 and IV.2, 


I f{x 


n \ 
NM) 


~ f* I — I/( x 7Vm) — InM^nm) + /wm(^jvm) — /iVM ~ f* + /nm\ 

* Tv + W' + (IV.7) 


where A^ fM (A) is an upper bound on the optimization error after n iterations of 
optimization algorithm A are applied to (Pnm)- In this chapter we consider several 
different expressions for A^ M (A) under various assumptions about the optimization 
algorithm, and hence also about (Pnm)- Because it appears difficult to quantify the 
rate of convergence of the total error, as in Royset and Pee (2011) we concentrate on 
the rate of convergence of its upper bound in (IV.7). From (IV.7), it is clear that the 
rate of convergence of that bound gives a guaranteed minimum rate of convergence 
of the total error. 

Based on (IV.7), it is evident that different choices of N, M , and n may 
result in different bounds on the total error. Let b G N be the computational budget 
available to apply n iterations of A to (Pnm)- To make it clear that the choice of 
N, M, and n would typically depend on b , we write 1V&, M&, and rib- We refer to 
{(rib, Nb, Mb)}^!, with nb,Nb,Mb G N for all b G N, as a discretization policy. A 
discretization policy specifies the level of discretization of (P) as well as the number 
of iterations of the optimization algorithm to complete for any computational budget. 
If rib, Nb, Mb —* oo, as b —» oo, then based on Assumption IV.1, the bound on the 
discretization error vanishes. If we assume that a convergent optimization algorithm 
is applied to solve (Pnm), then the optimization error and, hence, the corresponding 
bound could vanish as well. For a given optimization algorithm A, and n, JV,MgN, 
we define the total error bound, denoted by e(A, n, N, M), as the right-hand side of 
(IV.7) and have 

e(A,n,N,M) 4 A + + A " KM (A). (IV.8) 
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In this chapter, we investigate the rate at which the total error bound e(A, n b , N b , M b ) 
vanishes as b tends to infinity for different discretization policies {{n b , N b , Mb)} b T 1 and 
optimization algorithms A. We give optimal discretization policies, which as discussed 
in detail below, attain the highest possible rate of convergence of the total error 
bound as the computational budget tends to infinity for a given class of optimization 
algorithms. 

The analysis that follows relies on the following assumption about the com¬ 
putational work needed by an optimization algorithm to complete n iterations on 
( Pnm )• 

Assumption IV. 3. There exist a a = a (A) G (0,oo), a p — fi(A) G (0, oo), and a 
v = v(A) G (0, oo) such that the computational work required by a given optimization 
algorithm A to complete n 6 N iterations on (Pnm), N, M E N, is no larger than 
cmN^M u . □ 

Assumption IV.3 holds with p = 1 and v — 2 if the optimization algorithm A is 
a steepest descent method. This is because each iteration of these algorithms requires 
the calculation of /nm(x) at the current iterate x E X N , as well as the evaluation 
of the gradient V /nm(x), which dominates the computational budget. In order to 
evaluate the gradient, using problem (GTP^ m ) from Chapter III as an example, 
we see from (III.60) that O(N) operations are required to find the physical states, 
while (III.61) and (III.73), respectively, indicate that 0(NM 2 ) operations are required 
to find the information state and adjoint. This leads to an overall computational 
complexity of 0(NM 2 ) to evaluate V /nm(x). 

The SQP algorithm in the TOMLAB SNOPT solver (see Gill et ah, 2007) also 
seems to empirically follow Assumption IV.3 with p = 1 and v — 2 when used to solve 
(i GTP C ) with K = 1 and L = 10 using Algorithm V.2, which will be discussed further 
in Chapter V, with p 0 = (7 t/4 , 0 G M^ 0 ), target and HVU parameters values given in 
Table 6 (except v rnax = 60), and searcher and detection rate parameters given in the 
k = 1 columns of Tables 7 and 8, respectively. Table 1 shows the computational time 
in seconds required to complete five iterations of the SQP algorithm for various values 
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of N and M. Table 1 also shows the fitted values for a zero intercept linear regression 
done in R, version 2.13.0 (see Team, 2010), for the linear model Y = aNM' 2 . Figure 
2 illustrates the close agreement between the actual and fitted values given in Table 
1 . 


Table 1. Actual and fitted computational time in seconds for five iterations of the 
SQP algorithm in the TOMLAB SNOPT solver. 


Instance 

N 

M 

Actual 

Fitted 

1 

16 

7 

1.854 

1.427 

2 

16 

9 

2.164 

2.359 

3 

32 

7 

2.620 

2.855 

4 

32 

9 

4.176 

4.719 

5 

64 

7 

4.532 

5.709 

6 

64 

9 

8.655 

9.438 

7 

128 

7 

10.152 

11.419 

8 

64 

11 

13.337 

14.099 

9 

128 

9 

15.279 

18.876 

10 

128 

11 

24.154 

28.197 

11 

128 

15 

42.168 

52.433 

12 

128 

17 

60.439 

67.348 

13 

256 

11 

69.144 

56.395 

14 

256 

15 

102.924 

104.866 

15 

256 

17 

120.854 

134.695 

16 

512 

15 

196.666 

209.733 

17 

512 

17 

289.129 

269.390 


In some implementations, particularly those involving discretization of opti¬ 
mal control problems using collocation methods, increasing the time discretization 
parameter, N, gives rise to a large number of equality constraints. This large number 
of equality constraints provides special structure that allows for sparse implementa¬ 
tions where the computational work may grow slowly as N increases. Assumption 
IV.3 may account for the slow growth of the computational work in N, by selecting 
a value for /i e (0,1). 

Based on Assumption IV.3, as in Royset and Pee (2011), we refer to a dis¬ 
cretization policy {(rib, A&, M b )} < ^ =1 as asymptotically admissible if aribN^M^/b —* 1, 
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-Actual 

-Fitted 


Figure 2. Regression fit for Y = aNM 2 model. 

as b —> oo. It is obvious that an asymptotically admissible discretization policy satis¬ 
fies the computational budget in the limit as b tends to infinity. In the next section, we 
determine optimal asymptotically admissible discretization policies and correspond¬ 
ing rates of convergence of the total error bound under different assumptions about 
the optimization algorithm and, therefore, the optimization error bound A 

C. RATE ANALYSIS FOR CLASSES OF ALGORITHMS 

From (IV.7) we see that the total error bound consists of discretization and 
optimization error bounds. The discretization error bound depends on the discretiza¬ 
tion levels N and M, but not on the optimization algorithm used. The optimization 
error bound depends on the rate of convergence of the optimization algorithm used 
to solve ( Pnm )• In this section, we consider four cases: First, we assume that the 
optimization algorithm solves (Pnm) in a finite number of iterations. Second, we 
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investigate optimization algorithms with a superlinear rate of convergence towards 
an optimal solution of ( Pnm )• Third, we consider linearly convergent optimization 
algorithms. Fourth, we deal with sublinearly convergent algorithms. 

1. Finite Optimization Algorithm 

Suppose that the optimization algorithm used to solve (Pnm) is guaranteed to 
find an optimal solution in a finite number of iterations, independently of N and M. 
We define finitely convergent algorithms based on the following definition adapted 
from Royset and Pee (2011). We note that our definitions for superlinearly, linearly, 
and sublinearly convergent algorithms are also adapted from Royset and Pee (2011). 

Definition IV. 1. An optimization algorithm A converges finitely on {(Pvm)})^^ m=m 
when Xjy M is nonempty for N > N, M > M, and there exists a constant n e N such 
that for all N > N, M > M, N, Me N, a sequence {xnm}™=o generated by A when 
applied to (Pnm) satisfies G X* NM for all n >n. □ 

No optimization algorithm converges finitely on {(Pnm)}™’^ m=m without 
strong structural assumptions on Xn and f nm(- ), such as linearity. In this disserta¬ 
tion, we are not concerned with instances of (Pnm) that are linear programs, which 
may allow for finite convergence. We include this case here to serve as an “ideal” case. 

It will be shown below that this case gives an upper bound on the rate of conver¬ 
gence of the total error bound using any optimization algorithm. Based on Definition 
IV.1, a finitely convergent optimization algorithm A fimte on {(Pnm)}^’^ m=m ^ ias 110 
optimization error after performing a large enough number of iterations. We define 
AWA fimte ) = 0 and e(„4. fimte , n, N, M) = Ah + Ah for n > n, N > N, and M > M, 
where K is as in Assumption IV. 1 and h, N, and M are as in Definition IV. 1. The 
next theorem gives the rate of convergence of the total error bound for this case. 

Theorem IV.4. Suppose that Assumption IV. 1 holds and that A^ nite is a finitely con¬ 
vergent algorithm on {(Pnm)}™’^ m= m> Wl ^ 1 ^ an> ^ ^ as Assumption IV. 1 and 
number of required iterations n as in Definition IV. 1. Suppose also that satis¬ 

fies Assumptions IV.2 and IV.3. If {(nt,, N b , M b )}'jfi =1 is an asymptotically admissible 
discretization policy with n b = ii for all be N, then 
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1 


lim inf log e(A finite , n bl Nb, M b ) > 


(IV.9) 


b^o o log b (h/P + v /<?) ’ 

where p and q are as in Assumption IV. 1, and p and v are as in Assumption IV.3. 
Furthermore, if N b /b 1 ^ tJ,+pu / q ' > — > cq G (0, oo) as b —>■ oo ? i/ien 

log e(A finite , n b , N b , M b ) 1 


lim inf 

6—>-oo 


(IV.10) 


log 6 (h/p + v/q)' 

Proof. For sufficiently large b, A^ M (^ fimte ) = 0 and e(v4 fimte , n, N, M) = jV + A-, 
where K is as in Assumption IV. 1. Consequently, for sufficiently large b, 


log e(A fimte , n b , N b , M b ) = log 

, ( f K K 

> log max {—=, —, 


(IV.ll) 


V 


NT Ml 


= max l log log 


K 


Ml 


(IV.12) 


Hence, 

log e(A fimte , n b , N b , M b ) 
log b 


> max 


-p log N b + log K —q log M b + log K 1 


log b 


log b 


I 


. (IV. 13) 


We now consider two cases, one for both of the terms inside the max function of 
(IV. 13) when it is greater than or equal to the other term inside the max function. 
Initially, we want the first term to be greater than or equal to the second term. 


This will be true if 




P log N b > —q log M b . 

(IV. 14) 

We note that (IV. 14) implies 




Ml > Nl u/q . 

(IV. 15) 
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If we then use the bound on obtained from (IV. 15), we have that 


log e(^I finite , n b , N b , M b ) 
log b 


> 


—p log N b + log K 
log b 

-p log N b + log K 


log 


lo S (vhvfAf ) + lo S 

-p log N b + log K 

lo s + log(ari) + log [ Nb 


pv+qpL 
Q 


-P 


log A- 
log N b 


crn 


(IV.16) 


log N b 


log (an) , pu+qp 
log N b "T q 


Next, we consider the case where the second term inside the max function of 
(IV. 13) is greater than or equal to the first term inside the max function. This will 
be true if 


-glog Ah, > -plog N b . 


(IV. 17) 


We note that (IV.17) implies 


K > M^ /p . 


(IV.18) 


If we then use the bound on N b obtained from (IV.18), we have 
log e(*4. fimte , n b , N b , M b ) -q log M b + log K 


log b 


> 


log b 

-q log M b + log K 


> 


lo s (vuw ) + log ( aflN b M b 

-q log M b + log K 


log 


crnNgM" 


+ log (cm) + log ( M t 


pv+qpL 


- Q + l °s K 

1 ^ log M b 


log 


(IV.19) 


lN b M b J I log(g-n) , pu+qp 
log M b log M b p 


log e(^4 fimte , n b , N b , M b ) 
log b 


We now consider 



as b —* oo. Since anNjfM%/b —> 1 as b —>■ oo, we consider two cases: one where at 
least one of the parameters N b or M b does not increase without bound as b —» oo, 
and another where iVj, —> oo and Mb —» oo as 6 —* oo. For the first case, at least 
one of the parameters iVj, or M h remains hnite as b —> oo. Then based on (IV.11), 
e(*4. fimte , n b , N b , M b ) G ^0, lim sup^^ , and we have that 


\oge(A {inite ,n b ,N b ,M b ) 1 

Inn inf --- - - > 0 > - 7 —- rT , 

b^o o lOgO (P/P+V/q) 


(IV.20) 


For the second case, where both of the parameters N b and M b increase without bound 
as b —> 00 , we have that 


lim - nf log e(A &nite ,n b ,N b ,M b ) 

b- 5>oo log 6 


> min < lim inf 

b—>00 


-n + lp g A ' 
n ^ log JVj, 


anN^M" 

O O 


log At, 


1 log(<m) , pu+qp 
log At, " r 9 


lim inf 

6 —>-oo 


= mm 


log X 
log M;, 


anN^Mt 
_ o o 


log M;, 

1 


1 log(tm) , pu+qp 
' log M b _l_ p 


{fi/p + v/qY (fi/p + v/q) 


{n/p + v/qY 


which completes the proof for the first part of the theorem. 

Next, let {(n b , N b , M b )}^ =1 be an asymptotically admissible discretization pol¬ 
icy with n b = n satisfying N b /b l ^ i - i+pv / q ' > —> a\ G (0, 00 ) as b —» 00 . For notational 
simplihcation we define 


B 


1 — 


1 

Kb v/v+v/q 

A ~p 


(IV.21) 
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Then, for sufficiently large b , e(A fimte , n, N, M) = and algebraic manipulation 

gives 


e(A &n * e ,n b ,N b ,M b ) 


b \<i/ v 


TS" / b \ 

= .Bi &* i / p+I '/ ? H - crnN^ ) 

' b \<l/ v Mi 


-l 

= b^i vJrV i q 


V crnN^ / 

( b Ka q / v n q / v N m / v 1 

B 1 + ,/ _ b^T q 


Mi 


bi/ h 


Algebraic manipulation can be used to show that 


(IV.22) 


, 1 n 

foliv+pv 2 /q l/^b = 6 U = 1. 


(IV.23) 


Then, since an b N{fM%/b —y 1, the sum of terms in brackets in (IV.22), with n, N, 
and M replaced by n, N b , and M b , respectively, tends to a constant as b —» oo. The 
conclusion for the second part of the theorem then follows from (IV.22) after taking 
logarithms, dividing by log b, and taking limits. □ 

Theorem IV.4 shows that e(^4 fimte , n b , N b , M b ) converges at a rate fr-VlAp+V 1 ?) 
under the stated discretization policy. From (IV.8) we see that the total error bound 
includes the discretization error bound. This means the total error bound cannot 
converge faster than the rate no matter which optimization algorithm is 

used. 

It would be difficult to implement the asymptotically admissible discretization 
policy stated in Theorem IV.4 because n may be unknown. Despite this difficulty, 
the rate of convergence obtained in Theorem IV.4 is useful as an upper bound on 
the rate that can be achieved by any optimization algorithm, and will be used as a 
benchmark for comparison in the case of superlinear, linear, and sublinear algorithms 
below. 

2. Superlinear Optimization Algorithm 

We will now consider superlinearly convergent optimization algorithms. We 
define superlinearly convergent algorithms based on the following definition. 
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Definition IV. 2. An optimization algorithm A converges superlinearly with order 
7 G (l,oo) on {(Pnm)}^’^ m= m w h en i s nonempty for N > N, M > M, and 

there exist constants n G N , c G [ 0 , oo), and p G [ 0 , 1 ) such that c 1 ^ 7_ 1 } {JnmFnm) ~ 
f*NM) < P and 


Jnm (^jym) Inm 
JnM ( X NM ) — /atm ] 7 


< C, 


(IV.24) 


for all n > n, n G N, N > N, and M > M, N, M G N. 


□ 


Definition IV.2 requires that the optimization algorithm achieve a superlinear 
rate of convergence for sufficiently large n. This is usually the case with 7 = 2 
for Newtonian methods applied to strongly convex instances of ( Pnm) with twice 
Lipschitz continuously differentiable objective functions. The next lemma gives a 
total error bound for a superlinearly convergent algorithm. 


Lemma IV.5. Suppose that Assumption IV. 1 holds and that yf su P er is a superlinearly 
convergent algorithm with orders G (1, 00 ) on {(Pnm)}^’^ m =m> ™’£/i N and M as 
in Assumption IV. 1. Let be the iterates generated by 7 f su P er when applied 

to (Pnm), NeN, N>N, Me N, M> M. Suppose also that 7 f su P er satisfies 
Assumptions IV.2 and IV.3. Then, there exist constants c G (0,1); K G [0, 00 ), and 
n G N such that 




n \ 
NMJ 


- F < A n n + 


I< 

Np + 


K 

Wp 


(IV.25) 


for all n > n, n G N, N > N, N G N, and M > M, M G N, where K, p, and q are 
as in Assumption IV. 1. 


Proof. The proof follows the same arguments as those for the proof of Lemma 1 in 
Royset & Pee (2011). □ 

From Lemma IV.5, we adopt the upper bound on the optimization error 


AW-4 supen A 


) — C K 


(IV.26) 


for a superlinearly convergent optimization algorithm 7 f su P er on {(Pnm)}^F^ m= m% 
where c and n are as in Lemma IV.5. Then for n, N, M G N, we define the total error 


bound 


e(^4 super , n, N, M) A d" K + A + A. 
v ’ ’ ’ ’ Np Mi 


(IV.27) 


101 



The next result shows that if we select a particular discretization policy, then a super- 
linearly convergent optimization algorithm results in the same rate of convergence of 
the total error bound as a finitely convergent algorithm. Therefore, the policy speci¬ 
fied in the following theorem is optimal in the sense that no other policy guarantees 
a better rate of convergence. 


Theorem IV.6. Suppose that Assumption IV. 1 holds and that *4. super is a superlin- 
early convergent algorithm of order 7 G (1, 00 ) on {(Pnm)}^’^ m= m> with N and 
M as in Assumption IV.1. Suppose also that *4. super satisfies Assumptions IV.2 and 
IV.3. If {(nb, N b , M b )}'jfL 1 is an asymptotically admissible discretization policy with 


or 


and 


then 


nb log7 v , , M . , 

= l and. -;-< 1, 


log ( Xi, 

n b log 7 


log 3 


log b 


- logo 


N h 


6 log 7 


rlog(d 


log b 


log c y 


pq + up 


—> Ox G (1, 00 ), 


—>■ 02 G (0, 00 ), 


pq+vp 


lim 

b—>00 


log e(^4 super , n b ,N b , M b ) 


(IV.28) 


(IV.29) 


(IV.30) 


, , - -n- 7 T) (IV.31) 

logo {V/P + v/q) 

where p and q are as in Assumption IV. 1, and p and u are as in Assumption IV.3. 


Proof. For notational simplification we define 


and 


n log ■7 


B^expIlogK + logel-^W^ 

- log c 


K 


Bo = 


b log 7 


rl°g(d 


log b 


log c , 


pq 

pq+vp 


Np 


B 3 


01 og 7 



-pq 

pq+vp 


(IV.32) 


(IV.33) 


(IV.34) 


102 



Algebraic manipulation of (IV.27) gives 


e(A super ,n b ,N b ,M b ) 

/ n log 7 

= exp I log k + ")/ los ( -iogc) 


log ( logt) ' 

10 &7V-logc, 


logc 


— B\ + -B 2 -B 3 + 


K 

(_b_\ q / v 
V cmN / 


b y/ u 

crnNP ) 

Mi 


+ -B 2 -B 3 + 


K 

( b y/ v 

\ anN^ ) 


(_b_y/ u 

V anNP ) 

Mi 


= b 3 


B 1 
£3 


+ £2 + 


crnNL 1 ) 9/ Ktj‘>l''nil v NV<il'' 


Mi 




£3 


(IV.35) 


We now focus on the following portion of the first term inside the square 
brackets in (IV.35), 

/ 

log b 


Bib^+vp = exp logc 


logc 


log 


n log -7 

log b 

- log c 


pq 

Qpq+vp 


(IV.36) 


To simplify the analysis, we take the logarithm and obtain 

pq 


log ft + log c 

= log ft + log b 
Since an b NjfM%/b —» 1, 


n log 7 

log b \ io g ( 


logc 


— log c ) 


+ 


log b 


pq + up 

n log 7 _ 


l0SC „..„ (logt) 1 "^)" + 


pq 


(— log c) log (-°°gc) 


pq + up 


(IV.37) 


n b log 7 



N b 


b log 7 

CTl °s(^f Sh) 


q 

pq+up 


—t Q 1 , 


—)■ 02 , 


and, due to the facts that logc < 0 and a\ — 1 and ^ < 1 or a\ e (1, 00 ), the 

expression in (IV.37) goes to — 00 , as b —» 00 , we see that the sum of terms in brackets 
in (IV.35), with n, N, and M replaced by n b , N b , and M b , respectively, tends to a 
constant as b —» 00 . The conclusion then follows from (IV.35) after taking logarithms, 
dividing by log b , and taking limits. □ 
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While no discretization policy can guarantee a rate of convergence better than 
l)-'/(jj/p+'-'/q)_ p j s possible to do worse. If. for example, {(rib, N b , M b )}^ =1 is an asymp¬ 
totically admissible discretization policy with 


w v , , pq ^, 

—» oi = 1 and - ; - < 1, 


b 1 / 2 


/iq + vp 


(IV.38) 


or 


and 


n b 

b 1 / 2 

N b 


b 1 / 2 


—y cl i G (1, oo), 


—1 ci '2 G (0, oo), 


pq+up 


(IV.39) 


(IV.40) 


then it can be shown using arguments similar to those used in the proof of Theorem 
IV.6 that 

log e(A super , n b , N b , M b ) 1 


lim 

6—>-oo 


log b 


2/x ■ 2v 
p q 


(IV.41) 


which is worse than the optimal rate by a factor of two. 


3. Linear Optimization Algorithm 

We now consider linearly convergent optimization algorithms. We define lin¬ 
early convergent algorithms based on the following definition. 


Definition IV. 3. An optimization algorithm A converges linearly on {( Pnmm=m 
when X* N is nonempty for N > N, M > M and there exist constants n G N and 
c G (0,1) such that 


Jnm ( x nai) /. 


NM 


< c, 


fNM ( X NM ) /iVM 

for all n > n, n G N, N > N, N G N, and M > M , M G N. 


(IV.42) 

□ 


Definition IV.3 slightly extends a standard definition of linear convergence to 
require that the rate of convergence coefficient holds for all N and M sufficiently large. 
This is satisfied by many gradient methods, such as the steepest descent method and 
projected gradient method, when applied to (Pnm) under the assumption that the 
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objective function is strongly convex and twice continuously differentiable and that 


X is convex. The next lemma gives a total error bound for a linearly convergent 
algorithm. 

Lemma IV. 7. Suppose that Assumption IV. 1 holds and that *4. lmear is a linearly 
convergent algorithm on {(Pnm)}™’^ m=m> w% lh N and M as in Assumption IV.1. 
Let {x n NM }™ =Q be the iterates generated by v 4 hnear when applied to ( Pnm ), N e N, 
N > N, M e N, M > M. Suppose also that there exists a constant CeR such that 
/nm^x 1 ^) < C for all n e N, N > N, N e N, and M > M, M e N, and that *4. lmear 
satisfies Assumption IV.2. Then, there exists a constant k e [0, oo) such that 


f(x 


n \ 
NMJ 


- f * < C n K, + 


K 

m + 


K 

Wv 


(IV.43) 


for all n > n, N > N, and M > M, where c and n are as in Definition IV. 3, and 
K, p, and q are as in Assumption IV. 1. 


Proof. Based on Assumption IV. 1 and the fact that v 4 linear is linearly convergent, we 
obtain that 


\f( X NAl) f I — Inm(x% m ) + jy p + fNM 


< c n ~ n 

< c" 


n * K K 

[./' j V A/f ( X jy j^f j fN m\ ~V ]yp J[pQ 


c 


K K \ \ K K 

O' — /* H-1 -11 H- — +-• (IV.44) 

J Np Mi Np Mi y ’ 


Hence, the results hold with k = c n (C — f* + -^ + □ 

From Lemma IV.7, we adopt the upper bound on the optimization error 


A A 


(IV.45) 


for a linearly convergent optimization algorithm *4. lmear on {(P/vm)}^^ m=m j w h ere 
c and k are as in Lemma IV.7. Then for n, N, M e N, we define the total error bound 

e( A M “’,n,N,M)± ffV + A + V (IV.46) 

The next result shows that if we select a particular discretization policy, then a 
linearly convergent optimization algorithm can also attain the best possible rate of 
convergence of the total error bound given in Theorems IV.4 and IV.6. 
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Theorem IV.8. Suppose that A hnear satisfies the assumptions of Lemma IV. 7 and, in 
addition, Assumption IV.3 holds. If {(n b , N b , M b )}£f i is an asymptotically admissible 
discretization policy with 


-n b log c 
log b 


—y q i £ 


PI 

pq + up : 



(IV.47) 


and 


N b 

q 

Hq+vp 


—6 logo 
cr log b 


—y 0,2 G (0, oo), 


(IV.48) 


where c is as in Definition IV.3, p and q are as in Assumption IV. 1, and p and v are 
as in Assumption IV. 3, then 


Hm log e( v 4. linear , rib, Nb, M b ) = 1 

oo log 6 p/p + v/q 


(IV.49) 


Proof. For notational simplification we define 


B i = 


K 


pq 

—blogc pq+vp 
<j log b 


Np 


and 


B 2 


/— felogc\ 
\ a log b ) 


-pq 

pq+vp 


Algebraic manipulation of (IV.46) gives 


(IV.50) 


(IV.51) 


Since 


e(A^,n b ,N h ,M b ) 

( —n log c log 6 ^ , D d , K (wvO 

exp log k H --—= logc + B x B 2 + 


b \q/ v 


log b — log c 


Tl log C 

exp I log k + -—— log 6 ) + BiB 2 + 


( b \i/ u Mi 

V crnNP ) 

K (^-) q/u 

\ anNv ) 


— n log c + B\B 2 
b log b 


log b J 1 z ( b \i/ v Mq 

\ crnNP ) 

K (->-) q/v 

V anN^ ) 


= Bo 


K 


— n log c 

B 2 b lo « 6 


B i 


( b \l/ v M q 

V onN^ ) 

( b y/ 1 ' j£ a q/v n q/v]ynq/v 


k cjtiNp , 


B 2 Mi 


bi/ v 


-n b log c 
log b 


—y Gp G 


pq 


pq + vp 


,oo), 


(IV.52) 
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the expression 


(IV.53) 


K 

— n log c 

B 2 b lo s 6 

goes to 0, as b —* oo. In addition, we know that crn b N(fM£/b —> 1, and 


N b 

Q 

Hq+vp 


— 61 ogc 
cr log b 


—> 0 , 2 , 


as 6 —> oo. Consequently, the sum of terms in brackets in (IV.52), with n, N, and 
M replaced by n b , N b , and M b , respectively, tends to a constant as b —> oo. The 
conclusion then follows from (IV.52) after taking logarithms, dividing by log&, and 
taking limits. □ 


4. Sublinear Rate of Convergence 

Next we consider sublinearly convergent optimization algorithms. We define 

sublinearly convergent algorithms based on the following definition. 

Definition IV.4. An optimization algorithm A converges sublinearly with degree 
7 G (0, oo) on {(Pnm)}^’^ m= m wh en X* N is nonempty for N > N, M > M and 
there exists a constant C G [0, oo) such that 

f nm (*£jvm) — f.X.U — (IV.54) 

for all n G N, N > N, N G N, and M > M, M G N. □ 


The subgradient method is sublinearly convergent in the sense of Definition 
IV.4 with 7 = 1/2 and C = D x Lf when ( Pnm ) is convex, where D x is the diameter of 
X N and Lf is a Lipschitz constant of / Am (•) on V; see Nesterov (2004), pp. 142-143. 
Based on Definition IV.4, we define the optimization error bound for a sublinearly 
convergent optimization algorithm yf sublm by 

^«(-4- ub,i ") & T, (IV.55) 

and the total error bound for n, N, M G N by 

a, N, M) 4 g + A + A. (IV.56) 
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The next theorem gives an optimal discretization policy for a sublinearly convergent 
optimization algorithm as well as its corresponding rate of convergence of the total 
error bound. 


Theorem IV.9. Suppose that Assumption IV. 1 holds and that A suhhn is a sublinearly 
convergent algorithm with degree 7 G (0, 00 ) on {( Pnmm=m’ w% ^ 1 N and M as 
in Assumption IV. 1. Suppose also that 7 f sublm satisfies Assumptions IV.2 and IV.3, 
and that {(rift, AT b , M b )} r fiL l is an asymptotically admissible discretization policy. Then, 


log e ( 7 f sublin , n b , N b , M b ) 
hm - - - - - 

b — ^00 log 0 


1 

I _ i _ H E ’ 

7 p q 


(IV.57) 


where 7 is as in Definition IV.fi p and q are as in Assumption IV. 1, and /j and v 
are as in Assumption IV.3. 

Furthermore, if n b /b 1 / ( ' tJ/ y/p+ l "y/ q + 1 ) —>• oq G (0, 00 ) and 


N b 


p.-y/p+v'y/q 
b P-l /P+v-y/ q+1 


pq+vp 


—t 0,2 G (0, 00 ), 


(IV.58) 


as b —» 00 , then 


lim 

b—>00 


log e(^4 sublm , n b , N b , M b ) 


log b 


I _|_ M _|_ v ' 


7 V 


Proof. For any n, N, M G N, 


C K 


K 


log e(^ sublin , n, N, M) = log — + — + — 
v y \rP Np Mi 

f (C K K 

“ I rV Np ’ M ( i 


K 


K 


= max <j log c - 7 log n , log —, log — 


Hence, 


log e(v4 sublin , n, N, M) 


> max 


log b 

-7 log n + log C —p log N + log K —q log M + log K 


(IV.59) 


(IV.60) 


(IV.61) 


log b ' log b ’ log b J ^ ^ 

We now consider three cases, one for each of the three terms inside the max function 
of (IV.62) when it is greater than or equal to the other two terms inside the max 
function. 
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Initially, we want the first term to be the greater than or equal to the other 
two terms. This will be true if 



—7 log n b > -p log N b 

(IV.63) 

and 




-7 log ri h > ~q log M b . 

(IV.64) 

We note that (IV.63) implies 




K > 

(IV.65) 

and (IV.64) implies 




> npr 

(IV.66) 


If we then use the bounds on N p and obtained from (IV.65) and (IV. 66 ), respec¬ 
tively, we have that 

log e(v4. sublm , n, N , M) —7 log n b + log C 

log b ~ log b 

-7 log n b + log C 

lo S ( vnJtMs ) + lo S ( an b N bK) 

> _ -7 log n b + log C _ 

r \ / pq+qitl+pv'y \ 

lo s {^kw) + loga + log p w ) 

_ 7 + MV 

= —---^^-. (IV.67) 

S \ !Trl b N b M b ) 1 logo- pq+qp-y+pu-y 

log n b ' log n b pq 

We now consider the case where the second term inside the max function of 
(IV.62) is greater than or equal to the other two terms. This will be true if 

—p log N b > ^7 log rift, (IV. 68 ) 

and 

— plogN b > —q\ogM b . (IV.69) 

We note that (IV. 68 ) implies 

n b > N P J\ (IV.70) 
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and (IV.69) implies 


Ml > Nl u/q . (IV. 71) 

If we then use the bounds on n b and Ml obtained from (IV.70) and (IV.71), respec¬ 
tively, we have that 

-p log N b + log K 


log e(.A sub , n, N, M) 
log b 


> 


log b 

-p log N b + log K 


> 


lo g { <rn b N^M^ ) + lo g K. 

-p log N b + log K 


lo g ( *n b N b mz ) + log^ + ^g (W{ 


pq+qpy+pis'-y 
<11 


~P 


log K 
log N b 


log 


™b N b M b 


(IV. 72) 


log Nb 


I logo- pq+q^+pv-l 

' log N h 97 


Finally, we consider the case where the third term inside the max function of 
(IV.62) is greater than or equal to the other two terms. This will be true if 



—q log M b > -7 log n b , 

(IV. 73) 

and 




—q log M b > p log N b . 

(IV. 74) 

We note that (IV.73) implies 




n b > M q J\ 

(IV. 75) 

and (IV. 74) implies 




Nl > M^ /p . 

(IV. 76) 
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If we then use the bounds on n b and N{f obtained from (IV.75) and (IV.76), respec¬ 
tively, we have that 


log e( v 4 sublin , n, N, M) 
log b 


> 


> 


—q log M b + log K 
log b 

-q log M b + log K 

lo S ( *n b N£Mz ) + lo S (<rn b N£MZ) 
-q log M b + log K 


lQ S +1 ° gff + l0g ( 


pq+qp^+pv^ 

M, 


-q + 


log K 
logM 6 


log 


log M b 


I logo- I 
" r log M b ^ 


pq+qp'y+pvy 

PI 


(IV. 77) 


We now consider 

log e(»A sublm , n bl N bl M b ) 
log b 

as b —)■ oo. Since an b N^M£/b — > 1, as b —> oo, we consider two cases: one where at 
least one of the parameters n b , N b , or M b does not increase without bound as b —> oo, 
and another where n b — > oo, N b —* oo and M b —* oo, as b —> oo. For the first case, at 
least one of the parameters n bl N b , or M b remains finite as b —> oo. Then based on 
(IV.60), e(A subhn ,n b ,N b ,M b ) G (o, limsup^^ yr + jr? + jjq \ and we have that 

\ o o V1 b J 


lim inf 

b—>oo 


log e(^4 sublin , n b , N b , M b ) 
log b 


> 0 > 


I I P _I_ V 
7 ' p g 


(IV. 78) 


For the second case, where all three of the parameters n b , N bl and M b increase without 
bound as b —* oo, we have that 

log e(*4 sublin , n b , N b , M b ) 


lim inf 

b—>oo 


log b 


> min < lim inf 

b—>oo 


-ry _J_ iggg 

' log n b 


™b N b M b J , logo- I pq+qp.'y+pv'y 

log n b log n b pq 
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logiVf, 
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, lim inf 

b—>oo 
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log A' 
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I logo- I pq+qm+PVl 
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7 p q ) 


1 P _I_ V 

7 p q 


(IV. 79) 


Which completes the proof for the first part of the theorem. 

Next, let {(n-6, W, Mfe)}^ =1 be an asymptotically admissible discretization pol¬ 
icy satisfying n b /b 1 / ( ' fJ "y/p+ l "y/ q + 1 ) —> e (0, oo) and 

N b 


It'l/p+V'y/q ^ pq+vp 
b M"Y/ P+v-Y / <J +1 


—y Cl 2 £ (0, oo), 


as b —>■ oo. For notational simplification we define 
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Then, algebraic manipulation of (IV.56) gives 

e(A suhhl \n b ,N b ,M b ) 

= B^^ 1 +B 2 B 3 + 
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b \i/ u 


\ crnNV ) 


= B, 


D ~7 

1 , M7 , n 
- h p q 
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+ b 2 + 


( b \i/ v M q 

\ crnN^ ) 

(dw) 9/l/ Ka^n^N^ 


B 3 M q 


b q / v 


(IV.80) 


(IV.81) 


(IV.82) 


(IV.83) 


(IV.84) 
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Algebraic manipulation can be used to show that 


MV | V'l 

b i^^ b v^ = b ° = i, 


(IV.85) 


MV + ^V o M7 | U 1 

g/v P ~ q _ pq P q _ pq 

b “ +V + 1 b “ + 7 + 1 b- q/v b “ + 7 + 1 m+vv = b 0 = 1. (IV. 86 ) 

Then, since an b N^M^/b — » 1, the sum of terms in brackets in (IV.84), with n, N, 
and M replaced by n^, V 6 , and M b , respectively, tends to a constant as b —» oo. The 
conclusion for the second part of the theorem then follows from (IV.84) after taking 
logarithms, dividing by log b, and taking limits. □ 

Based on Theorem IV.9, the rate of convergence of the total error bound for a 
sublinearly convergent optimization algorithm is worse than that which is achievable 
by finite, superlinear, and linear algorithms (see Theorems IV.4, IV. 6 , and IV. 8 ), 
even when the optimal discretization policy given in the second part of the theorem 
is used. We also note that if 7 tends to infinity, then the rate for the sublinear case, 
when the optimal discretization policy is used, tends to that of the finite, superlinear, 
and linear cases as it should. 


D. CONCLUSIONS 

In this chapter, we consider the rate of convergence of a bound on the error 
between the objective function evaluated at iterates generated from the discretized 
problems and the optimal value of the original problem as a computational budget 
b tends to infinity. We see that in the case of superlinear and linear optimization 
algorithms, the best possible rate of convergence is ^/(m/p+V?^ where /i and v are 
positive parameters related to the computational work per iteration in the optimiza¬ 
tion algorithms, and p and q are positive parameters related to the error in the 
numerical methods used to construct the finite-dimensional problems. We identify 
specific optimal discretization policies for both the superlinear and linear cases that 
achieve this best possible rate. It is impossible to improve upon this rate due to 
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the fact that there is always some level of discretization error. If other discretization 
policies are utilized, it is possible that they will result in substantially slower rates of 
convergence. For sublinear optimization algorithms, with optimization error bounded 
by C/ri 7 , C > 0 , 7 > 0 , after n iterations, the best possible rate of convergence is 
&- 1 /( 1 / 7 +/Vp+"/g)_ This rate is slower than what can be achieved with either a super- 
linear or linear optimization algorithm. The optimal discretization policy specified 
in Section IV.C.4 achieves the rate 5 ~ 1 /( 1 /t'+Vp+ ! V | 2 ). Table 2 gives a summary of the 
results for the different optimization algorithms we consider in this chapter. 


Table 2. Comparison for optimization algorithms. 


Optimization Algorithm 

Asymptotic rate of decay of error bound 

Finite 

-1 

b^/p+v/q 

Superlinear 

-1 

b^/p+v/q ( w ith optimal policy) 

Linear 

=T 

Ijp/p+v/q ( w ith optimal policy) 

Sublinear 

-1 

bi/'i+p/p+v/q (with optimal policy) 


The results from this chapter provide insight regarding the type of numeri¬ 
cal method used to solve the differential equations when implementing a discretiza¬ 
tion algorithm to solve one of the generalized optimal control problems discussed in 
Chapter III. If a linear or super linear optimization algorithm A is used to solve the 
finite-dimensional optimal control problems, with Euler’s method used to numerically 
approximate the solution of the differential equations and Simpson’s rule used to nu¬ 
merically approximate the spatial integration, then p — 1 and q = 4 and we have 
e(A,n, N, M) = A^ M (A) + ^ + jji- If we assume that Assumption IV.3 holds for 
the reasons discussed in Section IV.B with p = 1 and v — 2, then the best possible 
rate of convergence is b ~ 2 / 3 . 

If a Runge-Kutta algorithm is used instead of Euler’s method to numerically 
solve the differential equations, it is clear from an analysis of Runge-Kutta algorithms 
(see, for example, Nagle & Saff, 1989 pp. 133-134) that although additional function 
and gradient evaluations are required, Assumption IV.3 could still hold with p = 1 
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and v — 2. If a second-order Runge-Kutta method is used instead of Euler’s method 
to numerically solve the differential equations, then p = 2 and q — 4 and we have 
e(A, n, N, M) = A r f KfM (A) + jpz + jfi- The resulting best possible rate of convergence 
is then b _1 . If a fourth-order Runge-Kutta method is used instead of Euler’s method 
to numerically solve the differential equations, then p — 4 and q — 4 and we have 
e(A,n,N,M) = A^ M (A) + jjs + jp ■ The best possible rate of convergence is then 
6 -4 / 3 . 

In order to establish an upper bound for the type of numerical method used 
to solve the differential equations, we let p —* oo. Then we have e(A,n, N, M) = 
A jy M ( A) + jj±. If we assume that Assumption IV.3 still holds with p = 1 and u — 2, 
then the best possible rate of convergence is b~ 2 . The analysis of this chapter therefore 
indicates that higher-order methods such as Runge-Kutta or pseudo-spectral may lead 
to better overall performance. However, this potential improvement in performance 
comes at the cost of additional complications in both the theory and implementation. 

We can also establish an upper bound for the type of numerical integration 
method used to evaluate the spatial integration if we let q —* oo. Then we have 
e(A,n, N, M) = A r ^ M (A) + If we assume that Assumption IV.3 still holds with 
p = 1 and v — 2, then the best possible rate of convergence is b~ p . The value of p 
depends on the choice of numerical method used to solve the differential equations, 
with p = 1 for Euler’s method, p — 2 for a second-order Runge-Kutta method, and 
p — 4 for a fourth-order Runge-Kutta method. Table 3 gives a summary of the 
results for the different numerical methods used to solve the differential equations 
and evaluate the spatial integration. 

The last row in Table 3 gives the asymptotic rate of decay of the error bound 
assuming “Ideal” methods are used to solve the differential equations as well as eval¬ 
uate the spatial integration. The term “Ideal” method means we let p —> oo and 
q —> oo. Then we have e(A, n, N, M) = A^ M (A). There is no benefit associated with 
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increasing N or M, so b — n and the resulting asymptotic rates are c 76 and c b for 
superlinear and linear optimization algorithms, respectively. 


Table 3. Comparison for numerical methods used to solve differential equations and 
evaluate the spatial integration. The optimization algorithm can be finitely, super- 
linear ly, or linearly convergent. The last row in the table gives the asymptotic rate 
of decay of the error bound assuming “Ideal” methods are used to solve the differ¬ 
ential equations as well as evaluate the spatial integration. The rates given are for a 
superlinear optimization algorithm with order 7 G (1, 00 ) and c G (0,1), and a linear 
optimization algorithm with rate constant c G (0,1). 


Numerical Method for 

Numerical Method for 

Asymptotic rate of 

Differential Equations 

Spatial Integration 

decay of error bound 

Euler 

Simpson 


2 nd Order Runge-Kutta 

Simpson 

b~ l 

4 th Order Runge-Kutta 

Simpson 

b=r 

“Ideal” p —> 00 

Simpson 

Ir 1 

Euler 

“Ideal” q —> 00 

b- 1 

2 nd Order Runge-Kutta 

“Ideal” q — > 00 

b~ 2 

4 th Order Runge-Kutta 

“Ideal” q — > 00 

b~ 4 

“Ideal” p —> 00 

“Ideal” q — > 00 

'yb — h 

c 1 or c 


We leave verification of the theoretical rate of convergence results developed 
in this chapter via numerical examples for future work. We note that this verification 
will be challenging clue to the fact that it will be difficult to determine what the 
optimal value of the objective function should be for the problems under considera¬ 
tion. The theoretical results of this chapter are still useful, however, as they provide 
insight regarding the choice of numerical methods for approximately solving the dif¬ 
ferential equations and approximating the spatial integrals as well as the optimization 
algorithms that can be utilized to implement a discretization algorithm that solves a 
generalized optimal control problem. 
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V. ALGORITHMS AND NUMERICAL 

RESULTS 

This chapter bridges the gap between the theory developed in Chapter III and 
implementable algorithms that can be used to solve the problems ( GTP ), ( GTP C ), 
(< GTP e ), ( GTP c ’ e ), (ITP e ), (ITP c ' e ), ( ITP p ), and ( ITP C ^). Again, for the sake 
of brevity we omit explicit treatment of the problems ( GTP e ) and ( GTP c,e ) in this 
chapter. The results developed in this chapter for ( ITP e ) and (JTP c,e ) can again 
be trivially extended to include ( GTP e ) and ( GTP c,e ), respectively, with a 1 replaced 
by a and (p l (•) replaced by 0(-). In Section V.A, we discuss problems defined on a 
real Euclidean space of coefficients needed to construct implementable algorithms, 
followed by a statement of the algorithms and their proofs of convergence, where 
appropriate. Then, in Section V.B.l, we use one of these algorithms, with a fixed 
discretization scheme, to solve an instance of ( GTP C ). This solution is then used to 
derive operational insights on how to better defend a HVU against an attack from a 
collection of small boat adversaries. We also give examples of numerical solutions to 
instances of ( ITP c,e ) and ( ITP C,P ) using fixed discretization schemes. In an effort to 
develop faster solution methods, in Section V.B.2 we use an algorithm with an adap¬ 
tive precision-adjustment scheme to solve ( GTP C ) and compare the results acheived to 
those obtained using different fixed discretization schemes. Finally, in Section V.B.3, 
we compare three heuristic methods designed for use onboard unmanned systems to 
provide solutions to ( GTP C ) in real time. 

A. IMPLEMENTABLE ALGORITHMS 

The approximating problems ( GTP NM ), (GTP^ M ), (ITP^ M ), (ITP^ e M ), 
(ITP^ m ), and (ITP^ p M ) in Chapter III were all defined using the function space H^. 
We would like to use existing nonlinear programming algorithms to solve these ap¬ 
proximating generalized optimal control problems. The issue with using these existing 
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algorithms is that they are defined only on real Euclidean spaces. For this reason, 
we need to define equivalent problem formulations in a corresponding real Euclidean 
space of coefficients that we can use to calculate numerical solutions. Equivalence is 
based on correspondence between solutions as in sections 4.3 and 5.6 of Polak (1997). 
We adopt the notation used in sections 4.3 and 5.6 of Polak (1997) to define these 
problems, and mirror the development found on pages 722-723 of Polak (1997) in this 
section. 

If we let e 3 denote the jth unit vector in M m , then the functions e J 7T J v A (-), 
j = 1 = 0, ...,N — 1, with Hn,k defined as in (111.54), form an orthonormal 

basis for the subspace Ln defined in (111.55). This means that Hn is in one-to-one 
correspondence with the real Euclidean space 

H n = R n x R mN . (V.l) 


The elements of H N are q = (£, u), with £ G M n and u = (h 0 , hq,..., un-i) G R m xl m x 
... x M m . Because Hn is in one-to-one correspondence with Hn, any q = (£, u) G Hn, 
with u(-) = T,^Wn,k{-)> corresponds to a unique fj = (£,h) G Hn, with u split 
into N, m-dimensional blocks u = (ho, hi, ...,un-i)- Then for any N G A7, we define 
the linear, invertible maps Wn : Hn —f Hn by 


N -1 


W N U-E 'U j k,'^N,k ,(*) J — (£5 ^1? •••? UN— l)) • 


K=0 


For any 77 , 77 ' G Hn and r] = Wn{j]) ) rf = Wn(v')i 

(v, v')h 2 = {£> O + [ («(*), “'(*)> dt 


N-l 


(£,£') + (^u k ,VNu', 


K=0 


N-l 


(£, O + U «) - fa h) Hn > 


Av—0 


(V.2) 


(V.3) 


where (•, •) denotes the Euclidean inner product on and (•, -)g denotes the Eu¬ 
clidean inner product on Hn- Which means for any q G Hn and fj = Wn i'q), 

INI h 2 = (v, v ) H2 = \\v\\s N = (v, v)h n • ( V - 4 ) 
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This shows that Wn (•) is a norm-preserving, or isometric, map whose adjoint is 
which means given any q G Hn and q G Hn, 

?))„„ = (try (V.5) 

Given an q = (£,u) G T/jv, if we let q = IV4 ( 77 ), then (III.60) defines the 

I _ >| v 

sequence of vectors 4 — 1 3% (j) \ i n for all A: = 1, 2 , K, by the recursion 

l J j =0 

(0 + 1)/JV) - if (j/N) = /h* (4‘0/JV). Vnu'u/n)) . 

i 6 {0,1.JV— 1} ,*5r(0) = e- (V. 6 ) 

Similarly, (III.61) and (III. 181) define the sequence of vectors zn — {zl{j]a)} N G K 

_ r_ 1 n 

and z l N = l 4 (j; a 1 ) [ Gl by the recursions 
l j j =0 

I< L 

z V n(U + 1 )/ N 'i a ) -z v N (j/ N ',a) = vEE^ (j/N),y l {j/N- a)), 

v fc=i ;=i 

J e {0,l,...,iV-l},4(0;a) = 0, (V.7) 


and 


z]f l (( 3 + i )/^;«') - 4 * ( i /^;«0 = ^ { x N k U/ N ),y l { j / N \ «*)) , 

v fc=i 


3 G {0,l,...,iV-l},4(0;^)=0, (V.8) 


respectively. 

Using the discretized “information state” given by the recursion (V.7), we 
define the approximate objective functions Jnm ■ H N —> M for any q G Hn and 
N G J\f by 

= hi (exp [- 4 (!; •)] 0 ( 0 ) > (V.9) 

and the corresponding discrete generalized optimal control problems by 

(GTP nm ) min f NM {fj), (V.10) 

v&h n 
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and 


(< GTP c nm ) min f NM (fj), (V.ll) 

Ci jv 

where 


Hc,w — 


x {u = (u 0 ,Ul; ...,U N -l) 


nmN 


Un G 


U,j G {0,1,..., IV — 1}} , (V.12) 


with U as in (III.9). 

Similarly, using the discretized “information state” given by the recursion 
(V.8), we define the approximate objective functions ip e NM '■ Hn —>• M and i/j P NM : 
H n —y M. for any fj G H N and N G AT by 

L 


^nm(v) = ( exp </> z (■))> 


(V.13) 


Z=1 


and 


nm(,v ) = n Im ( exp ■ 


Z=1 


(V. 14) 


We dehne the corresponding discrete generalized optimal control problems by 


(ITP l nm ) 

min 

rj£H N 


(V. 15) 

(ITP% e M ) 

min 


(V.16) 

(ITP p nm ) 

min 

rj£H N 

^Pnm(j} S ) i 

(V.17) 

(ITP^ n p m ) 

min 

>teH CjAr 

i’NM (v)- 

(V.18) 

Based on Exercise 4.3.7 in Polak 

(1997), the problems (GTPnm, 

). (gtp c nm ), 


(. ITP e NM ), (ITP C j^ m ), ( ITP p nm ), and {ITP C ^ M ) are equivalent to the problems 
(GTP nm ), (CrTP c NM ), (), (/TP£ e M ), and respectively, 
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in the sense that they are related by a nonsingular linear transformation, and that at 
corresponding points the values of their corresponding optimality functions are the 
same. 

We now state generic algorithm models based on Master Algorithm Model 
3.3.12 from Polak (1997) that can be used to solve the problems ( GTP ), (GTP C ), 
(. ITP e ), ( ITP c,e ), ( ITP P ), and ( ITP C,P ). Our algorithms define the “outer” it¬ 
erations, while the “inner” iterations are defined by user supplied iteration maps 
A((P),fj), where A denotes one iteration of a nonlinear programming solver applied 
to problem (P) starting from rj. We begin with the unconstrained problems and state 
the following algorithm, which is based on a fixed discretization scheme. 

Algorithm V.l. Approximately solves (GTP). 

Data: N 0 G Af, M 0 G N 3 x M 3 ' 5 , r / 0 G H^ o . 

Step 0. Set N = N 0 , M = M 0 , and fjo = Wn(vo)- 

Step 1. Generate {fji}™ 0 using fj i+1 G A((GTP NM ),fji). □ 

Algorithm V.l also approximately solves ( ITP e ) and (ITP P ), with (GTP^m) re- 
placed by ( ITP e NM ) and (ITP p nm ), respectively. For the constrained problems, we 
have the following algorithm. 

Algorithm V.2. Approximately solves (GTP C ). 

Data: N 0 G A/", M 0 G N3 x M3, r ] 0 G H,, jA r 0 . 

Step 0. Set N = N 0 , M = M 0 , and fjo = fFv(r/o)- 

Step 1. Generate {ry,:},A 0 using fj i+1 G A((GTP c NM ),fji). □ 

Algorithm V.2 also approximately solves (ITP c ' e ) and (ITP C,P ), with (GTP C NM ) re¬ 
placed by (ITP C ^ M ) and (ITP C ]^ M ), respectively. 

In order to improve the run-time performance of Algorithms V.l and V.2, we 
also develop algorithms based on an adaptive precision-adjustment scheme. Algo¬ 
rithms based on adaptive precision-adjustment keep the discretization parameters, 
N and M, small when the algorithm starts and is “far” from the optimal solution. 

5 Recall that N 3 = {m € 2N + 1| m > 3}, as defined in Proposition III.26. 
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We use the optimality function to determine how “far” the current iterate is from 
the optimal solution. The values of N and M are then increased according to suc¬ 
cessor functions as the algorithm gets closer to the optimal solution. One advantage 
these algorithms have over fixed precision algorithms is that the computational work 
remains relatively low until the algorithm gets close to the optimal solution. An¬ 
other advantage is that solutions obtained using low levels of discretization are used 
to “warm start” the algorithm at higher discretization levels. Algorithms based on 
adaptive precision-adjustment fit nicely within the framework of consistent approxi¬ 
mations because we know that as the values of N and M tend to infinity, the solution 
of the approximating problem converges to the solution of the original problem. These 
algorithms use successor function k : Af —>■ Af defined by 

k(N) E {N' E Af\N' > N} . (V.19) 

It is clear from (V.19) that k(N) is an integer larger than N e Af. We begin with 
the unconstrained problems and state the following algorithm, which is based on an 
adaptive precision-adjustment scheme. 

Algorithm V.3. Solves ( GTP ). 

Data: N 0 e Af, M(N) :A/’->N 3 xN 3 , r) 0 e H^. 

Parameter: f3> 0. 

Step 0. Set i — 0, N — N 0 , M = M(N 0 ), and f), = Wjv(^i)- 
Step 1 . Compute an fj i+1 G A((GTP NM ( N )),fji). 

Step 2. If 9(W^(fj i+1 )) > —(Pi/N + /3 2 /(Mi(N)) 4 + ^ 3 /(M 2 (A‘)) 4 ), set V * N = 
W^ 1 , and replace N by n(N). 

Step 3. Replace i by i + 1, and go to Step 1. □ 

Algorithm V.3 also solves ( ITP e ) and ( ITP P ), with (GTP nm(n)) replaced by 

(ITP e NM ( N ^) and ( ITP p nm , n A , and d(-) replaced by 6 e (-) and 0 P (-), respectively. Then 

the following theorem is a direct consequence of Theorem III.29. 

Theorem V.l. If the sequences and {rj^} are constructed by Algorithm V.3, 

then 
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(i) if the sequence {rf N } is finite, then the sequence { 77 *}°^ has no accumulation 
points; and 

(ii) if the sequence {r]* N } is infinite, then every accumulation point rf G H° of 

{v*n}^ o’ satis fi es 9(v*) = 0. □ 

When Algorithm V.3 is used to solve ( ITP e ) and ( ITP P ), a result similar to Theorem 
V.l with Of) replaced by 6 e f) and 0 p f), respectively, follows as a direct consequence 
of Theorem III.52. 

For the constrained problems, we have the following algorithm. 

Algorithm V.4. Solves ( GTP C ). 

Data: Nq G A/", M(N) : J\f —> M 3 x M 3,770 £ H c jv 0 . 

Parameter: /A,/3 2 , fis > 0. 

Step 0. Set i — 0, N — N 0 , M = M(N 0 ), and fy = Wjy(r]i). 

Step 1. Compute an g i+ 1 G A((GTP c NM ^),fji). 

Step 2. If 0 c (W^(fj i+l )) > —(fii/N + (3 2 /(Mi(N)) 4 + /3 3 /(M 2 (iV)) 4 ), set rj* N = 
Wfif 1 (fji+i) , and replace N by k,(N). 

Step 3. Replace i by i + 1, and go to Step 1. □ 

Algorithm V.4 also solves (JTP c,e ) and ( ITP C ’ P ), with ( GTP° NM replaced by 
{ITP c f^ M ^ N fj and {ITP c f^ M ^ N f), and 6 c f ) replaced by 0 c,e f) and 0 c,p f), respectively. 
Then the following theorem is a direct consequence of Theorem III.30. 

Theorem V.2. If the sequences {fji}°l 0 and { 77 ^} are constructed by Algorithm V.f, 
then 

(i) if the sequence { 77 ^} is finite, then the sequence {fii}fiL 0 has no accumulation 
points; and 

(ii) if the sequence { 77 ^} is infinite, then every accumulation point rj* G H c of 

{'77tv} i=o’ satisfies 9 c {rf) = 0. □ 

When Algorithm V.4 is used to solve ( ITP c,e ) and ( ITP C,P ), a result similar to The¬ 
orem V.2 with 0 c f) replaced by 0 c,e (-) and 0 c,p f), respectively, follows as a direct 
consequence of Theorem III.53. 
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Theorems V.l and V.2 are the end result of the theory we develop in Chapter 
III. Theorems V.l and V.2 establish the fact that the discretization schemes we 
develop in Chapter III result in implementable algorithms that are guaranteed to 
converge to stationary points of the original problems. 

B. NUMERICAL RESULTS 

1. Fixed Discretization Schemes 

In this section we present numerical results based on the situational description 
given in Chapter II. We start by solving three problem instances given in Table 4 
using Algorithm V.2, with parameter values given in Table 5. The column headings 
K and L in Table 4 denote the number of searchers and targets, respectively. The 
column heading fjo in Table 5 represents the initial “guess” for the algorithm. The 
first three elements of fjo are the initial headings of the searchers measured from the 
horizontal axis. The initial control portion of fj 0 is 0, which represents the zero vector 
of length Nq. We note that we do not include the initial position of the searchers in 
fjo, as they are assumed to be fixed and not part of the decision vector. The initial 
positions for the searchers are provided in Table 7. The parameter values given in 
the first row of Table 5 are used with Algorithm V.2 for ProbA and ProbC, and the 
parameter values given in the second row are used with Algorithm V.2 for ProbB. 

Table 4. Fixed discretization problem instances. 


Instance 

Problem Class 

K 

L 

ProbA 

(< GTP C ) 

3 

10 

ProbB 

(. ITP c ’ e ) 

3 

2 

ProbC 

(ITP C,P ) 

3 

1 


Table 5. Algorithm V.2 parameters. 


No 

M 0 

Vo 

200 

(25,25) 

(tt/4, tt/2, tt/2, 0) 

320 

(25,25) 

(tt/4, tt/2, tt/2, 0) 


124 





In ProbA, ProbB, and ProbC we consider the case of a HVU operating in a 
two-dimensional area of interest (AOI), measuring 70 nautical miles (nm) by 70 nm. 
The HVU follows a straight line trajectory, with constant heading £3 and constant 
speed u° from its initial position x°(0). The parameter values for the HVU are given 
in Table 6 . While in transit, the HVU is under threat of attack from L targets. The 
target trajectories are conditioned upon the random vector, a, with realization in 
A C M 2 . For our numerical results, the target trajectories are conditioned upon a 
random starting time between zero and one hour, and a random starting location on 
either of the vertical sides of the AOI. The numerical solutions we find for ProbA, 
ProbB, and ProbC are all based on a receding horizon search, where we are planning 
for the next hour. 

To generate numerical solutions for ProbA, ProbB, and ProbC, we implement 
and run Algorithm V.2 in MATLAB 7.11.0 (version 2010b) (see Math Works, 2011) 
on a 3.46 GHz PC with two quad-core processors, using Windows 7 Pro 64-bit, with 
24 GB of RAM. We use the SQP algorithm in TOMLAB SNOPT solver (see Gill 
et ah, 2007), with the default major and minor optimality tolerance as the stopping 
criteria. 

In ProbA, the HVU is under threat of attack from L = 10 targets. We assume 
that the distributions for the targets’ starting time and starting location are uni¬ 
form and independent. The range of the uniform distribution for starting location is 
[0,140], as we combine both sides of the AOI into a single segment. The range of the 
uniform distribution for the starting time is [0,1]. The targets follow deterministic 
trajectories, {y\t;a) : 0 < t < l}, given a, where y l (t]a ) = (y[(t; a ), y l 2 (t; a)) T G M 2 
is the position of the I th target at time t. 

The target trajectories are generated by solving the optimal control problem 
described in Section II.C with an additional constraint for each target. The additional 
constraint requires the targets to hit the HVU at a 90 degree angle of incidence at the 
final time, tf. The parameter values used to formulate this optimal control problem 
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are given in Table 6. The starting position for the I th target, y l 0 (a), is on one of the 
vertical sides of the AOI, and is determined based on the level of spatial discretization. 
We assume that the separation in starting position between members of the swarm is 
1.5 nm. The target trajectories are truncated because they are generated by solving 

Table 6. Target and HVU parameter values. The target parameter values are the 
same for all targets l — 1,2 ,..., L. 


v l o 

v f 

V‘min 

^max 

ylftar* 

x°(0) 

T .o 

v° 

15 kts 

19 kts 

1 kt 

35 kts 

250 (35,0) (nm,nm) 

tt/2 

25 kts 


a time-optimal control problem, and our search problem uses a receding horizon 
approach based on the planning horizon [0,1]. An example of this truncation is 
illustrated in Figure 3, where the trajectory of the HVU is indicated in blue, the first 
component of the vector M gives the number of starting locations for the swarm, 
and the second component of the vector M gives the number of starting times for 
the swarm. The graph on the top left of Figure 3 shows the trajectories for a single 
swarm consisting of L = 10 targets. The middle graph in the top row of Figure 3 
indicates the different trajectories that the swarm would follow if it started from five 
specific locations. The trajectories for the different starting locations are shown in 
alternating colors, red and then cyan. The graph on the top right of Figure 3 shows 
the different trajectories that the swarm would follow if it started from five specific 
locations and two specific starting times. The trajectories for the different starting 
locations are shown in alternating colors, red and then cyan, for the first starting time. 
The trajectories for the different starting locations are shown in alternating colors, 
green and then purple, for the second starting time. The graphs on the bottom row of 
Figure 3 are the same as their counterparts in the top row, except they are truncated 
to show how far the targets (and HVU) can progress on their desired trajectories 
during the 1 hour planning horizon. 

In an effort to detect the targets, the HVU has an escort consisting of K = 3 
searchers. The goal of the searchers is to minimize the probability of failing to detect 
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Figure 3. Full and truncated target trajectories for L = 10. Top Row Left: M = (1,1) 
Full, Top Row Middle: M = (5,1) Full, Top Row Right: M = (5, 2) Full. Bottom 
Row Left: M = (1,1) Truncated, Bottom Row Middle: M = (5,1) Truncated, 
Bottom Row Right: M = (5, 2) Truncated. 

any of the targets within the planning horizon [0,1]. For any t G [0,1], let x k {t ) = 
(x k (t), x k (t), x k (t)) T G M 3 be the physical state of the k th searcher at time t, where 
x\(t) G M and x k (t) G M are the horizontal and vertical components of the location 
of the k th searcher, respectively, and x k (t) G K is the heading of the k th searcher 
measured from the horizontal axis. We assume that the searchers are Dubins vehicles 
traveling at a constant velocity v k > 0. The control input for the k th searcher, u k G M, 
is the rate of change of the heading (also known as the yaw rate) of the k th searcher. 
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Then, the searcher dynamics described in general terms by (11.22) are given by 


h k (x k (t),u k (t)) 4 


v k cos x k (t) 
v k sin x k (t) 
u k (t) 


\ 


, k = 1,2,3. 


(V.20) 


We consider two different search platforms for our numerical results. One is modeled 
after an SH-60 helicopter, and the other is modeled after a guided missile destroyer 
(DDG). The searcher parameters are given in Table 7, where k — 1 is the helicopter 
and k — 2, 3 are the destroyers. The initial heading of the k th searcher is given by £ k , 
and the initial position of the k th searcher is given by x fc (0). The yaw rate limit for 
the DDG is based on discussions with CDR David Schiffman, USN, a surface warfare 
officer in the Operations Research Department at NPS (personal communication, 
November 2010), and the yaw rate limit for the helicopter searcher is based on the 
limit for the SH-60 helicopter found in A1-H60BB-NFM-000 (SH-60 NATOPS Flight 
Manual, 1996). 


Tabic 7. Searcher parameter values. 



k = 1 

k = 2 

k = 3 

e 

7r/4 

vr/2 

tt/2 

V k 

120 kts 

25 kts 

25 kts 

rr fc (0) 

(35,0) (nm,nm) 

(25,0) (nm,nm) 

(45,0) (nm,nm) 

Yaw Rate Limit 

1885 

hour 

250 A———■ 

hour 

250 ^ 

hour 


We let r k,l (x k (t),y l (t',a)) > 0 denote the detection rate at time t for the k th 
searcher at (x k (t),x k (t)) T when the I th target is located at y l (t ; a). For our numerical 
results, we adopt the Poisson Scan Model (see, e.g., p. 3-1 in Washburn, 2002) and 
set 


r^((^(t),4(t)) T ,^(t;a)) = A fc $[(F fc -p((^(t),4W) T ,l/ / (t;«)))/a fc ], (V.21) 


where $(•) is the standard normal cumulative distribution function, X k is the scan 
opportunity rate, F k is a sensor parameter, a k reflects the variability in the received 
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signal strength, and p k ((x k (t), x k (t)) T , y l (t; a)) is used to model the signal loss, which 
depends on the distance between the target and the searcher; see for example Figure 
4.5 on page 93 in Wagner et ah (1999). The detection rate functions for the helicopter 
searcher, shown in Figure 4, and the DDG searchers, shown in Figure 5, are based 
on p k ((x k (t),x k (t)) T ,y l (t;a)) = a k \\(x k (t), x k (t)) T — y l (t]a)\\ 2 + b k , with parameter 
values given in Table 8. The parameter values are chosen to reflect reasonable sensor 
performance based on discussions with other naval officers in the Operations Research 
Department at NPS who have served on these platforms (CDR David Schiffman, USN, 
and CDR Douglas Burton, USN, personal communication, November 2010). 


Table 8. Detection rate parameter values. 



k = 1 

k = 2 

k = 3 

x k 

1.1 

1 

1 

ppk 

90 

90 

90 

a k 

0.3 

0.3 

0.3 

b k 

20 

60 

60 

a k 

100 

15 

15 


The dimension of the decision vector, fj(-), in ProbA is KNq+K = 603, because 
we determine a control input for each of the searchers at every time step as well as the 
optimized initial heading for each searcher. Before we give our numerical results, we 
first show baseline results based on the current concept of operations (CONOPS) in 
use by the U.S. Navy. This CONOPS has the helicopter escort orbiting the HVU at a 
distance of 2.5 nm and the DDG escorts keeping station on parallel courses to the HVU 
at a distance of 10 nm. The CONOPS is based on discussions with helicopter pilots 
and surface warfare officers in the Operations Research Department at NPS (LCDR 
Harrison Schramm, USN, LCDR Ron Cappellini, USN, and CDR David Schiffman, 
USN, personal communication, February 2011). The resulting trajectories are shown 
in Figure 6. In Figure 6, the target trajectories are plotted in alternating colors of 
red and cyan to separate the specific starting positions of the swarm. It should be 
noted that only the set of target trajectories based on a starting time of t = 0 and all 
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Figure 4. Detection Rate Function for Helo Searcher. 

starting locations is plotted in Figure 6. The additional target trajectories based on 
other starting times between [0,1] are not shown in the interest of keeping the plot 
legible. The trajectory for the HVU is shown in blue, the DDG escort trajectories are 
shown in black, and the helicopter escort trajectory is shown in green. The objective 
value, which represents the probability that all of the searchers fail to detect any of 
the targets in [0,1], for this set of trajectories is 0.9589. It should be noted that as the 
targets get closer to the HVU the probability that all of the searchers fail to detect 
any of the targets will be lower, but this improvement might come too late for the 
HVU to effectively defend itself. 

The resulting trajectories from Algorithm V.2 are shown in Figure 7. The 
colors for the trajectories are the same as those described above for Figure 6. The 
helicopter trajectory is particularly good at illustrating how the algorithm gets the 
searchers to areas where they can accumulate as much probability “mass” as possible, 
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Figure 5. Detection Rate Function for DDG Searcher. 

and keeps them there to the extent allowed by the constant speed and turn rate 
constraints. The objective value is now 0.3954, which is significantly better than the 
baseline. This indicates there is a benefit to leaving the vicinity of the HVU and 
seeking out the potential attackers. The time required to solve ProbA was 26.82 
hours, including 1.54 hours to generate the target trajectories. 

In ProbB the HVU is under threat of attack from L — 2 targets. The target 
l = 1 is assumed to leave from the right-hand side of the area of interest and target 
l — 2 is assumed to leave from the left-hand side of the area of interest. The target 
trajectories are again conditioned upon a random starting time between zero and 
one hour, and a random starting location along the appropriate vertical side of the 
AOI. We assume that the target’s starting time and location are independent random 
variables. We assume that the distribution for starting position for both targets is 
uniform with range [0,70]. We assume that the distribution for starting time for 
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Figure 6. Trajectories based on current CONOPS. HVU trajectory is blue. Helicopter 
trajectory is green. DDG trajectories are black. Target trajectories alternate red and 
cyan. 

target l — 1 is triangular, with the target being twice as likely to leave at t = 0 as 
it is to leave at t — 1. We assume that the distribution for starting time for target 
l — 2 is uniform with range [0,1]. 

The target trajectories are generated in the same manner as they were for 
ProbA with the parameter values given in Table 9. 

Table 9. Target and HVU parameter values. The target parameter values are the 
same for targets l — 1,2. 


Vo 

V f 

V L min 

V-max 

LL 

x°(0) 

coo 

v° 

10 kts 

12 kts 

1 kt 

35 kts 

250 (35,0) (nm,nm) 

tt /2 

25 kts 


The HVU again has an escort consisting of K = 3 searchers. The goal of the 
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Figure 7. Trajectories based on Algorithm V.2 on ProbA. HVU trajectory is blue. He¬ 
licopter trajectory is green. DDG trajectories are black. Target trajectories alternate 
red and cyan. 

searchers is to maximize the expected number of targets detected within the planning 
horizon [0,1]. The searcher and detection rate parameters are the same as those given 
in Tables 7 and 8, respectively. 

The dimension of the decision vector, fj(-), in ProbB is again given by KN 0 + 
K = 963, which is different value than it was for ProbA due to the different parameter 
values used in Algorithm V.2. The resulting trajectories from Algorithm V.2 are 
shown in Figure 8. The color codes for the trajectories are the same as described 
above for Figure 6, except that target / = 1 trajectories are shown in red and target 
l = 2 trajectories are shown in magenta. Again, we note that only the set of target 
trajectories based on a starting time of t — 0 and all starting locations is plotted 
in Figure 8. The objective value, which represents the expected number of targets 
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detected in [0,1], for this set of trajectories is 0.4542. The time required to solve 
ProbB was 28.71 hours, including 0.33 hours to generate the target trajectories. 



Figure 8. Trajectories based on Algorithm V.2 on ProbB. HVU trajectory is blue. 
Helicopter trajectory is green. DDG trajectories are black. Target trajectories are 
magenta and red. 

In ProbC the HVU is under threat of attack from L — 1 target. The target is 
assumed to leave from one of the two vertical sides of the area of interest. The target 
trajectories are conditioned in the same way as they were for ProbA and ProbB. We 
assume that the distributions for starting time and starting location are independent 
and uniform. The range of the uniform distribution for starting location is [0,140], 
as we again combine both sides of the AOI into a single segment. The range of the 
uniform distribution for starting time is [0,1]. 

The target trajectories are generated in the same manner as they were for 
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ProbA with the parameter values given in Table 10. We note that we now consider a 
high-speed target, with a maximum speed of 60 kts. 


Table 10. Target and HVU parameter values. 


^0 

v f 

^min 

v 1 

max 

yl,td r 

rr 0 (0) 

, r o 

x 3 

v° 

10 kts 

12 kts 

1 kt 

60 kts 

250 (35,0) (nm,nm) 

vr/2 

25 kts 


The HVU again has an escort consisting of K = 3 searchers. The goal of the 
searcher is to minimize the probability of failing to detect any of the targets during 
the planning horizon [0,1]. The searcher and detection rate parameters are the same 
as those given in Tables 7 and 8, respectively. 

The dimension of the decision vector, fj(-), in ProbC is KN 0 + K = 603. The 
resulting trajectories from Algorithm V.2 are shown in Figure 9. The color codes for 
the trajectories are the same as described above for Figure 6, except that the target 
trajectories are shown in red. Again, we note that only the set of target trajectories 
based on a starting time of t = 0 and all starting locations is plotted in Figure 9. 
The objective value, which represents the probability that all of the searchers fail to 
detect any of the targets during [0,1], for this set of trajectories is 0.7979. The time 
required to solve ProbC was 2.59 hours, including 0.15 hours to generate the target 
trajectories. 

The results in this section show that Algorithm V.2 is tractable for as many 
as three searchers and ten targets. The results also indicate that it is beneficial 
for the searchers to leave the vicinity of the HVU and seek out potential threats in 
order to improve the overall probability of detection or expected number of targets 
detected. The results of ProbC illustrate the effect of a high-speed target. The 
helicopter searcher does not spend too much time on one side of the AOI, but rather 
attempts to cover what it can on one side and then flies over to the other side of 
the AOI. Finally, we see that obtaining numerical solutions for ProbA, ProbB, and 
ProbC using a fixed discretization scheme requires a significant amount of computing 
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Figure 9. Trajectories based on Algorithm V.2 on ProbC. HVU trajectory is blue. 
Helicopter trajectory is green. DDG trajectories are black. Target trajectories are 
red. 

time. There is evidence that using an adaptive discretization scheme can produce 
solutions that are equivalent to those obtained using fixed discretization schemes, 
but at a lower computational cost; see, for example, He & Polak (1990) and Section 
3.3.3 in Polak (1997). We consider an algorithm based on an adaptive discretization 
scheme in the next section. 

2. Adaptive Discretization Schemes 

In this section we obtain numerical solutions for ProbA using an adaptive 
discretization scheme based on Algorithm V.4 and compare them to results obtained 
for ProbA using the fixed discretization scheme of Algorithm V.2. The parameters 
/3i, @ 2 , and fa used to define what it means for the optimality function to be “close 
enough” to zero, the successor function, and the function M(N) used to specify 
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the level of spatial discretization as a function of the time discretization parameter 
considered in this section are summarized in Table 11. The column heading Set in 
Table 11 refers to the set of algorithm parameters used to obtain the solution for 
ProbA using either Algorithm V.2 or V.4, as indicated in the table. We note that 
the expressions for M t (N) given for Algorithm V.4 indicate that the initial spatial 
discretization parameters are (5, 5) and that they are each increased by adding 2 (for 
Sets 4 and 6) or 4 (for Sets 5 and 7) to the previous value every time the value of the 
time discretization parameter is doubled according to n{N). 


Table 11. Algorithm parameters used to obtain solutions for ProbA, where K = 3, 
L = 10, and 7 = 1,2. For all Sets, fj 0 = (vr/4, ir/2, ir/2, 0). 


Set 

Algo. 

N 0 

M 0 

Pi 

p2 

/i.3 

k(N) 

Mi(N) 

1 

V.2 

80 

(13,13) 

N/A 

N/A 

N/A 

N/A 

N/A 

2 

V.2 

120 

(17,17) 

N/A 

N/A 

N/A 

N/A 

N/A 

3 

V.2 

320 

(25,25) 

N/A 

N/A 

N/A 

N/A 

N/A 

4 

V.4 

10 

N/A 

10" 9 

5 x 10~ 10 

5 x 10~ 10 

N' = 2 N 

5 + 21og 2 ()))) 

5 

V.4 

10 

N/A 

10~ 9 

5 x 10~ 10 

5 x 10" 10 

N' = 2 N 

5 + 41og 2 (^) 

6 

V.4 

10 

N/A 

10" 8 

5 x 10” 9 

5 x 10” 9 

N' = 2 N 

5 + 21og 2 ()))) 

7 

V.4 

10 

N/A 

10~ 8 

5 x 10" 9 

5 x 10~ 9 

N' = 2 N 

5 + 41og 2 ()))) 


The solutions we obtain are based on the same scenario described in Section 
V.B.l for ProbA. The target trajectories are conditioned and generated as described 
above for ProbA using parameter values given in Table 6. The goal of the K = 3 
searchers is to minimize the probability that they all fail to detect any of the targets 
during [0,1]. The searcher and detection rate parameters are the same as those given 
in Tables 7 and 8, respectively. 

The data collected by using Algorithms V.2 and V.4 is shown in Figure 10. 
The objective function values for all seven algorithm parameter sets were evaluated 
using a “high precision 1 ' discretization level of N = 2560 and M = (91,91) to ensure 
they were all compared against a common standard, ft should be noted that when 
using the adaptive discretization scheme it is possible for the objective function value 
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to go up when the algorithm increases the values of the discretization parameters. An 
example of this behavior is seen for Set 6 after approximately 10000 seconds, where the 
algorithm increases N from 40 to 80 and M from (9,9) to (11,11). Based on Figure 
10, it is clear that all of the adaptive precision schemes eventually achieve a solution 
quality within 1.42% of the best solution obtained using one of the fixed discretization 
schemes. While the fixed discretization scheme for Set 2 does better than the adaptive 
discretization schemes after approximately 4000 seconds, it is important to note that 
all of the adaptive discretization schemes reach an objective value between 0.385 and 
0.39 more than a thousand seconds before the fixed discretization scheme for Set 2 
does. It is also important to note that all of the adaptive schemes perform better than 
the fixed schemes based on Sets 1 and 3. This implies that it is possible to do better 
for a specific problem if you happen to select a “good” fixed discretization scheme at 
the outset, but doing so without any prior knowledge about the problem is difficult. 
The adaptive schemes are fairly robust in the sense that they do well regardless of the 
initial choice of parameter values. The results in Figure 10 indicate that the adaptive 
precision schemes for Sets 5 and 7 offered the best overall performance compared to 
the other adaptive discretization schemes. 

Figure 11 shows the resulting trajectories from Algorithms V.2 and V.4. The 
color codes for the trajectories are the same as described above for Figure 6. Again, 
we note that only the set of target trajectories based on a starting time of t — 0 
and all starting locations is plotted in Figure 11. It is clear from the plots in Figure 
11 that the overall “shape” of the final searcher trajectories is the same regardless 
of the algorithm and parameters used to obtain the solution, with the caveat that 
the helicopter searcher trajectory is the mirror image when it goes to the opposite 
side of the AOI. The times used to plot the trajectories in Figure 11 for the fixed 
discretization schemes (Sets 1-3) were selected so that enough computational time 
had elapsed to allow the solution to stabilize. For the adaptive discretization schemes 
(Sets 4-7), the trajectories were only captured when the algorithm increased the 
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Figure 10. Comparison of fixed and adaptive precision schemes using Algorithms V.2 
and V.4, respectively. For Sets 1 and 2 computation was terminated between 10000 
and 15000 seconds because the solution had stabilized. Because Set 3 did not begin 
to stabilize until after 20000 seconds, the horizontal axis was extended. 


discretization level. Because of this limitation, the times used to plot the trajectories 
in Figure 11 for the adaptive discretization schemes were selected based on the largest 
computational time data that was available. 
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Figure 11. Trajectories for ProbA based on Algorithms V.2 and V.4. Top Row Left: 
Set 1 after 10339 seconds, N = 80, M = (13,13). Top Row Right: Set 2 after 
12201 seconds, N = 120, M = (17,17). Second Row Left: Set 3 after 25147 seconds, 
N = 320, M = (25,25). Second Row Right: Set 4 after 29852 seconds, N = 80, 
M = (11,11). Third Row Left: Set 5 after 23152 seconds, N = 40, M = (13,13). 
Third Row Right: Set 6 after 22167 seconds, N = 80, M = (11,11). Bottom Row 
Middle: Set 7 after 39998 seconds, N = 80, M = (17,17). 
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3. Real-Time Methods 

Although the adaptive discretization schemes discussed in Section V.B.2 offer 
improved run-time performance over fixed discretization schemes, they are still not 
fast enough to be implemented as real-time solution methods for use onboard a UAS 
or USV. In this section we consider three heuristic approaches that have the potential 
to be used as real-time solution methods onboard unmanned vehicles. This section 
begins with a discussion of how we implement the three heuristic methods. We then 
give a description of how the problem instances to be solved by these three methods 
were generated. Next we provide numerical results for the three methods. Finally, we 
state our conclusions based on the numerical results achieved by the three methods. 

a. Heuristic Methods 

We consider three different heuristic methods in this section that can 
be used to solve generalized optimal control problems. For simplicity we assume that 
we have only a single searcher, and we therefore omit the superscript notation for 
the searcher number. It should be noted that the development that follows could be 
trivially extended to include multiple searchers. The first method we consider uses 
Algorithm V.2. The second and third methods are based on fitting polynomials to 
determine optimal searcher trajectories. The first polynomial-based method attempts 
to optimize directly over the coefficients of the polynomials, while the second poly¬ 
nomial based method uses an indirect method similar to that found in Yakimenko 
(2000) and Ghabcheloo et al. (2009). The first polynomial-based method solves the 
following modified version of ( GTP C ). 


(i GTP 


c,dir\ . 


min 

a 


Mi M2 

/L W iJ eX P 

*=1 3 = 1 



* P( a i,j ) 
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s.t. V Q) 

o' 

W! 

;N- 1, j — 1,2 

(v. 22 ) 

■ (A\ 

x i XT 

2 \N / 

d i-1 

= (]y) , 7 = 0,..., N- 1, j = 1,2 

1=1 

(V.23) 

v — 0.05v 

< 

( ±l (A’ i2 (A) 

< v + 0.05u, 7 = 0 ,. 

.., A — 1 (V.24) 

a 

_ r 1 2 d. 12 

[a l5 CL^ • • • 5 Qj2' ) "" 


(V.25) 

°i 

= an( 0 ) 


(V.26) 

a 2 

= 3 : 2 ( 0 ), 


(V.27) 


where W tJ are the weights for Simpson’s rule, x (*) = ? Oii j are 

the discretization points at which the integrand is evaluated, d is the order of the 
polynomials, v is the searcher’s speed, and xi(0) and £ 2 ( 0 ) are the initial position of 
the searcher. We note that (V.24) constrains the searcher’s speed to be within 5% of 
v. 

We then use the following algorithm to solve (GTP c,dir ). 

Algorithm V.5. Approximately solves (GTP c,dir ). 

Data: N 0 e A f,M 0 eN 3 x N 3 , a 0 e M 2 x M 2 x ... x M 2 . 

Parameter: d G N, d > 1. 

Step 0. Set N = N 0 , M = M 0 , and a = ao . 

Step 1. Generate {a*}?2 0 using a i+ 1 G A((GT P c,dir ), a,). □ 

The second polynomial based method uses an indirect method to deter¬ 
mine the coefficients of the polynomials. We use the term indirect because instead of 
directly optimizing over the polynomial coefficients, we instead optimize over termi¬ 
nal conditions which can be used to obtain the polynomial coefficients. The method 
used is similar to that described in Yakimenko (2000) and Ghabcheloo et al. (2009), 
except that we omit deconffiction of trajectories. In our problem we also assume the 
searcher travels at constant velocity, so the searcher’s acceleration is zero over the 
entire planning horizon. The description of the method that follows is similar to that 
found in Yakimenko (2000) and Ghabcheloo et al. (2009), but is included here for 
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the sake of completeness. We denote by x(t) = {x\ (r) , x 2 (t)) t the desired path 
followed by the searcher, parameterized by the virtual arc r € [0 ,77], where 77 is the 
total virtual arc length between the initial and final positions of the searcher. We 
represent the state as a function of time using the same notation as the state as a 
function of the virtual arc length. The meaning should be clear from the context. We 
represent the coordinates x\ and X 2 by algebraic polynomials of degree d given by 

d 

x j (r) = J2 a J r ^ i = 1 ’ 2 > (V-28) 

i =0 

where r* indicates that r is raised to the i th power. The degree d of the polynomials 
Xj(t), j = 1 , 2 is related to the number of boundary conditions that must be satisfied. 
In the formulation that follows, we use the prime sign / to indicate d/dr, // to indicate 
the second derivative operator, and ill to indicate the third derivative operator. As 
before, we use dot symbols above variables to indicate derivatives with respect to 
time, t. Given 77 and terminal constraints Xj( 0 ), Xj( 0 ), x"( 0 ), x'"( 0 ), 37(77), Xj(t/), 
Xj (t y), and x'"(Tf), for j = 1,2, the coefficients of the seventh-order polynomial can 
be computed from 


(l 

0 

0 

0 

0 

0 

0 

0 ^ 


b?) 


* Zj(0) ^ 

0 

1 

0 

0 

0 

0 

0 

0 


a ) 


Xj( 0) 

0 

0 

2 

0 

0 

0 

0 

0 


a ) 


x?(0) 

0 

0 

0 

6 

0 

0 

0 

0 


a J 


x?(0) 

1 

T f 

T f 

r] 

T f 


-/ 6 

T f 


a 4 j 


x j( T f) 

0 

1 

2 Tf 

3 t) 

4 T f 

f 

6t / 



S 5 


x j( T f) 

0 

0 

2 

6 Tf 

12 r] 

20r3 

30 rf 

42 r} 


a J 


x "( x f) 

v° 

0 

0 

6 

24t/ 

60 rj 

120 r) 

210 rf) 


\ a V 


\ x 7( r f)/ 


(V.29) 

In order to express the relationship between r and t, we define A(r) = Then the 
temporal and spatial derivatives of x satisfy 

x = Xx'. (V.30) 
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As in Yakimenko (2000) and Ghabcheloo et al. (2009) we choose A to be an affine 
function of r, that is, A(r) = Ao + A; A ° r, with Ao = ||i(0)||, and A/ = ||i(l)||. 
By integrating t = A (r), the virtual arc r and time t are related by the following 
equations 


T f 


i°g A/Mo) ’ 


r 

r / 






\ f -\ 0 




if A/ = A 0 
if A/ ^ A 0 

if A/ = A 0 
if A/ 7 ^ A 0 


(V.31) 


(V.32) 


The second polynomial based method solves the following modified version of ( GTP C ), 
where Table 12 gives a summary of the relationship between the boundary conditions 
needed to solve (V.29) and the decision vector b. 


Table 12. Relationship between decision vector and boundary conditions for indirect 
polynomial method. 


b\ 

b 2 

h 

6 4 


h 

b 7 

h 


x z( T f) 

A (y) 

AAA 

<( 0 ) 


<(Tf) 

X 2(jf) 


jyc,indir\ . 

( Mi M 2 

min EE Wij exp 

b { mi j i 


N 


N -1 L 

EE 

7=0 1=1 


d 

s.t. Xj (t) = a i r *’ 3 = !, 2, Vr e [0, t/] 

i=0 

n — 0.05n < A(r) ||x'(r) || <v + 0.05n, Vr e [0, tj\ 

In problem (G'TP c,indM '), Wjj are the weights for Simpson’s rule, cqj are the discretiza¬ 
tion points at which the integrand is evaluated, d is the order of the polynomials, and 
v is the searcher’s speed. The total arc length r/ is computed using (V.31) and the 
spatial paths x{r) and speed profiles x are given by (V.29) and (V.30), respectively. 
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In order to discretize r, we set t — 7 = 0,1,... -,.N — 1 in (V.32). We note that 

for j = 1 , 2 , Xj( 0 ) are the starting coordinates for the searcher, a;"( 0 ) = x"(r/) = 0 , 
and x'j( 0 ) are determined by solving the initialization problem described in Algorithm 
V. 6 . 

We then use the following algorithm to solve (GTp c ^ indir y 
Algorithm V.6. Approximately solves (GT P c,indir ). 

Data: N 0 G A f , M 0 G N 3 x N 3 , rjo G H C; jv 0 - 

Parameter: d G N, d > 1. 

Step 0. Set N = N 0 , M = Ado, and fjo = IUv(r/o)- 

Step 1. Generate {? 7 i }™ =0 using rj i+ i G A((GTP c NM ), rii), where n is the first iterate 
satisfying a given optimality tolerance. 

Step 2 . Compute x l from fj n . 

Step 3. Set b 3 = xj( 1), j = 1 , 2 , b 3 = cos(a; 3 (l)), b 4 = sin(x 3 (l)), and b 3 = 0, 
3 = 5, 6 ,7, 8 . 

Step 4. Generate {h}^ using b l+ \ G A((GTP c ' mdir ), bi). □ 

In Step 1 of Algorithm V .6 we are using Algorithm V.2 to solve 
(■ GTP C NM ), and stopping it after satisfying a given optimality tolerance. In Step 
2 of Algorithm V .6 we use the solution, fj n , from Step 1 to compute the searcher’s 
trajectory. In Step 3 of Algorithm V .6 we set bj = Xj( 1), j = 1, 2, where (xi(l), 2 : 2 ( 1 )) 
is the final position of the searcher we calculated in Step 2 . We also use fj n from Step 2 
to determine the searcher’s final heading, which we resolve into components to obtain 
values for b 3 and 64. We initialize the remaining components of b to zero and move 
on to Step 4, where we determine an optimized value for b. 

b. Problem Instances 

All of the problem instances are based on a HVU operating in a 70 nm 
by 70 nm area of interest using a receding horizon approach, and a one hour planning 
horizon. For these instances, the HVU is under threat of attack from L = 10 targets. 
The target trajectories are conditioned and generated as described above for ProbA. 
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The parameter values used to generate the target trajectories are given in Table 6. 
We again assume that the separation in starting position between members of the 
swarm is 1.5 nm. We assume that the HVU has a single helicopter escort, whose 
goal is to minimize the probability that it fails to detect any of the targets during 
the planning horizon [0,1]. With the exception of Xi(0) and x 2 (0), which are given 
in Table 13, the searcher and detection rate parameters are the same as those given 
in Tables 7 and 8, respectively, with K — 1. 

We randomly generate the starting position for the searcher, given in 
columns two and three of Table 13. The xi(0) value was selected from the range 
[20,50], and the £ 2 ( 0 ) value was selected from the range [0,30]. We let the vertical 
starting position for the targets on the side (or sides) of the area of interest be a 
continuum between [0,140] and again randomly select the lower and upper bounds 
for the starting position of the targets. The generated values are given in columns 
four and five of Table 13. 

Table 13. Real-time problem instances. For all instances, problem class is ( GTP C ), 
K = 1, and L = 10. 


Instance 

2 h(0) 

£ 2 ( 0 ) 

Target LB 

Target UB 

Side Scenario 

ProbD 

34.56 

24.01 

19.86 

59.05 

LHS 

ProbE 

47.47 

11.77 

23.97 

91.77 

Both 

ProbF 

23.81 

27.4 

88.53 

134.05 

RHS 


c. Numerical Results 

In this section, we provide numerical results obtained using Algorithms 
V.2, V.5, and V.6 to solve the problem instances given in Table 13. We implement and 
run all Algorithms in MATLAB 7.11.0 on the same PC described in Section V.B.l, 
again using SNOPT with default major and minor optimality tolerance as the stop¬ 
ping criteria except as noted below. We separate the results for each problem instance 
into two tables, one based on the fixed discretization method given in Algorithm V.2, 
and one based on the polynomial methods given in Algorithms V.5 and V.6. We also 
provide a figure for each problem instance that shows the resulting trajectories from 
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the three solution methods. The result from Algorithm V.2 in Figures 12, 13, and 
14 uses the first set of parameters that achieve an objective function value smaller 
than either of the objective function values achieved by Algorithms V.5 and V. 6 . We 
plot all results using a value of IV = 320 and M = (25, 25) to ensure a common 
standard. This need to “level the playing held” creates two different issues based 
on the solution methodologies. The results given in Figures 12, 13, and 14 are from 
solutions generated with N values smaller than 320. As a result, the control inputs 
for the searcher resulting from Algorithm V.2 are linearly interpolated to N = 320 to 
generate the plots. Algorithms V.5 and V .6 generate polynomial coefficients, so there 
is no issue with the N = 320 used for plotting being larger than the N value used to 
generate the coefficients. There is an issue, however, when it comes to ensuring that 
the searcher maintains a constant speed equal to v along its trajectory. We generate 
the plots for Algorithms V.5 and V .6 using the path for the searcher defined by the 
polynomial coefficients, but we ensure that the distance traveled between each of the 
N = 320 searcher waypoints is limited to where v is the velocity of the searcher. 

Tables 14, 16, and 18 provide results from Algorithm V.2. The initial 
control input for the searcher is the zero vector of appropriate length, and the initial 
heading for the searcher is chosen based on the side scenario. For the left hand side 
and both side scenarios, the initial heading for the searcher is 37 t/4 . For the right 
hand side scenario, the initial heading for the searcher is 7 t/4 . Tables 14, 16, and 18 
list the values of N and M used to run Algorithm V.2, but all of the objective function 
values in row three of these tables are evaluated using a common discretization level of 
N = 320 and M = (25, 25). Row four in Tables 14, 16, and 18 gives the time required 
for the optimization only, and does not include time to build the target trajectories. 

To implement Algorithm V.5, we use d = 7. The initial input for the 
polynomial coefficients is chosen based on the side scenario. For the left hand side 
scenario, all of the coefficients are zero with the exception of a) which is set to a value 
of — aq(0). This results in an initial horizontal line trajectory for the searcher from its 
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starting location of (xi(0), x 2 (0)) and ending at (0, 24.01). For the both side scenario, 
we set a\ = —88.98, af = —0.37, a\ = 0.03, al = 45.64, a\ = —0.40, a\ = 0.04, 
and all the other coefficients equal to zero. This results in an initial straight line 
trajectory for the searcher with slope -0.507, beginning at (oq(0),a; 2 (0)) and ending 
at (—41.85,57.05). For the right hand side scenario, all of the coefficients are zero 
with the exception of a\ which is set to a value of 70 — aq(0) = 46.19. This results 
in an initial horizontal line trajectory for the searcher from its starting location of 
(aq(0),x 2 (0)) and ending at (70,27.4). Tables 15, 17, and 19 list the values of N 
and M used to run Algorithm V.5, but all of the objective function values in row 
four of these tables are evaluated using a common discretization level of N = 320 
and M = (25,25). Row five in Tables 15, 17, and 19 gives the time required for the 
optimization only, and does not include time to build the target trajectories. We note 
that for ProbD the stopping criteria for Algorithm V.5 is not the default major and 
minor optimality criteria for SNOPT. Because we are interested in real-time methods, 
we also set a maximum number of major iterations. In the case of ProbD, SNOPT 
reaches the 700 major iteration limit prior to satisfying its default optimality criteria. 
This is the reason why the searcher trajectory shown in middle plot of Figure 12 goes 
so far to the left, beyond the target trajectories. The algorithm was stopped before 
it could finish optimizing the trajectory to include as much probability “mass” as 
possible, while still satisfying the velocity constraints. 

To implement Algorithm V.6, we use d — 7. Because Algorithm V.6 is 
effectively initialized by Algorithm V.2, we again use the zero vector of appropriate 
length as the initial control input for the searcher. For the left hand side and both 
side scenarios, the initial heading for the searcher is 37 t/ 4. For the right hand side 
scenario, the initial heading for the searcher is 7 t/ 4. Tables 15, 17, and 19 list the 
values of N and M used to run Algorithm V.6, but all of the objective function 
values in row four of these tables are evaluated using a common discretization level of 
N = 320 and M = (25, 25). Row five in Tables 15, 17, and 19 gives the time required 
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for the initialization and optimization only, and does not include time to build the 
target trajectories. We note that the stopping criteria for Algorithm V.6 is not the 
default major and minor optimality criteria for SNOPT. Again, based on our desire 
to find real-time methods we set the stopping criteria as 15 major iterations. This is 
the reason why the searcher trajectories shown in the right plots of Figures 12 and 14 
go so far beyond the target trajectories. The algorithm was stopped before it could 
finish optimizing the trajectories to include as much probability “mass” as possible, 
while still satisfying the velocity constraints. We note that clue to the complicated 
relationship between the objective function and the decision vector, we use finite 
differences in SNOPT to estimate the gradients. Finally, we note that we initialize 
ProbE using N = 20 and M = (9,9) and not N = 10 and M = (5,5). This is 
because this method is highly dependent on a good initial input for the searcher’s 
final position, and the lower levels of discretization do not provide a solution where 
the searcher goes to both sides of the AOI to look for the targets. 
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Table 14. Algorithm V.2 results for ProbD. 


N 

10 

15 

20 

25 

30 

35 

40 

320 

M 

(5,5) 

(7,7) 

(9,9) 

(ii,ii) 

(13,13) 

(15,15) 

(15,15) 

(25,25) 

Obj. 

0.0476 

0.0294 

0.0324 

0.0206 

0.0097 

0.0118 

0.0092 

0.0088 

Time 

(sec.) 

10.67 

48.76 

141.77 

203.48 

361.53 

612.03 

884.68 

16223.48 


Table 15. Algorithm V.5 and V.6 results for ProbD. 



Direct 

Indirect 

N 

10 

10 

M 

(5,5) 

(5,5) 

Obj. 

0.0246 

0.0927 

Time (sec.) 

317.40 

176.50 



Figure 12. Trajectories for ProbD based on Algorithms V.2, V.5, and V.6. Left: 
After 203 seconds, Algorithm V.2, N = 25, M = (11,11) for solution, N = 320, M = 
(25,25) for plot. Middle: After 317 seconds, Algorithm V.5, N = 320, M = (25,25). 
Right: After 176 seconds, Algorithm V.6, N = 320, M = (25,25). 
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Table 16. Algorithm V.2 results for ProbE. 


N 

10 

15 

20 

25 

30 

35 

40 

320 

M 

(5,5) 

(7,7) 

(9,9) 

(ii,ii) 

(13,13) 

(15,15) 

(15,15) 

(25,25) 

Obj. 

0.5161 

0.5022 

0.2767 

0.2715 

0.2715 

0.2622 

0.2593 

0.2555 

Time 

(sec.) 

8.29 

32.60 

52.51 

93.55 

175.11 

284.54 

431.95 

43000.29 


Table 17. Algorithm V.5 and V.6 results for ProbE. 



Direct 

Indirect 

N 

10 

20 

M 

(5,5) 

(9,9) 

Obj. 

0.3526 

0.2918 

Time (sec.) 

66.05 

750.00 



*1 




Figure 13. Trajectories for ProbE based on Algorithms V.2, V.5, and V.6. Left: After 
53 seconds, Algorithm V.2, N = 20, M = (9,9) for solution, N = 320, M = (25, 25) 
for plot. Middle: After 66 seconds, Algorithm V.5, N = 320, M = (25,25). Right: 
After 750 seconds, Algorithm V.6, N = 320, M = (25,25). 
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Table 18. Algorithm V.2 results for ProbF. 


N 

10 

15 

20 

25 

30 

35 

40 

320 

M 

(5,5) 

(7,7) 

(9,9) 

(ii.ii) 

(13,13) 

(15,15) 

(15,15) 

(25,25) 

Obj. 

0.0553 

0.0330 

0.0246 

0.0265 

0.0249 

0.0237 

0.0234 

0.0216 

Time 

(sec.) 

9.29 

39.40 

74.77 

225.61 

464.91 

818.13 

1008.03 

35608.20 


Table 19. Algorithm V.5 and V.6 results for ProbF. 



Direct 

Indirect 

N 

10 

10 

M 

(5,5) 

(5,5) 

Obj. 

0.0734 

0.0421 

Time (sec.) 

67.85 

185.19 



Figure 14. Trajectories for ProbF based on Algorithms V.2, V.5, and V.6. Left: After 
39 seconds, Algorithm V.2, N — 15, M — (7, 7) for solution, N = 320, M = (25, 25) 
for plot. Middle: After 68 seconds, Algorithm V.5, N = 320, M = (25,25). Right: 
After 185 seconds, Algorithm V.6, N = 320, M = (25,25). 
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d. Conclusions 

Based on the results presented in this section, it appears that the 
method which is best suited to provide real-time solutions is Algorithm V.2. The 
fixed-precision method implemented in Algorithm V.2 consistently provides solutions 
with the best objective values and the lowest computational times. Table 20 summa¬ 
rizes the comparison between the different methods. To compute the values in Table 
20, we assume that the objective value obtained from Algorithm V.2 with N = 320 
and M = (25, 25) is the “correct” value. For ProbD, we use the objective values 
from Algorithm V.2 for N — 25 and M = (11,11) and N — 20 and M = (9,9) to 
compare with Algorithms V.5 and V.6, respectively. For ProbE, we use the objective 
value from Algorithm V.2 for N = 20 and M = (9, 9) to compare with Algorithms 
V.5 and V.6. For ProbF, we use the objective value from Algorithm V.2 for N = 15 
and M = (7, 7) to compare with Algorithms V.5 and V.6. The At columns in Table 
20 give the additional time in seconds necessary for the Algorithms V.5 and V.6 to 
compute their respective solutions. In fairness, it should again be stressed that the 
implementation of Algorithm V.6 utilized finite differences to estimate the required 
gradients. This puts Algorithm V.6 at a disadvantage over the others, although we 
attempt to mitigate this effect by stopping the algorithm after only 15 major itera¬ 
tions. Future studies should be done with the gradients computed explicitly to ensure 
a more balanced comparison. 


Table 20. Comparison of real-time methods. 


Algo. 

ProbD 
% Error 

At 

ProbE 
% Error 

At 

ProbF 
% Error 

At 

V.2 

134.8%/269.3% 

N/A 

8.29% 

N/A 

53.10% 

N/A 

V.5 

180.90% 

113.92 

37.98% 

13.55 

239.99% 

28.46 

V.6 

957.20% 

34.73 

14.19% 

697.49 

94.92% 

145.8 
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VI. 


CONCLUSIONS AND FUTURE WORK 


A. CONCLUSIONS 

Continuous timc-and-space search problems arise in numerous applications 
such as anti-submarine warfare, search and rescue operations, and protection of 
HVU’s from small boat attack. We consider the problem of detecting targets that 
seek to harm a HVU, and formulate it as a generalized optimal control problem. 
While the theory that we develop allows for more general cases, we assume that 
the targets follow deterministic trajectories, given realizations of random variables 
that provide information about their initial conditions. We consider two classes of 
problems, ( GTP ) and (/TP), based on whether the targets act in a coordinated or 
independent manner, respectively. For both classes, we consider objective functions 
that minimize the probability that all of the searchers fail to detect any of the targets 
during the planning horizon and maximize the expected number of targets detected 
during the planning horizon. 

We develop discretization schemes to solve (GTP) and ( ITP ) problems, and 
show that the finite dimensional problems are consistent approximations to their in¬ 
finite dimensional counterparts. We then use these discretization schemes to develop 
implementable algorithms that can be used to solve the problems (GTP), ( GTP C ), 
(< GTP e ), (GTP c,e ), (ITP P ), ( ITP C ' P ), (ITP e ), and (ITP c,e ). We provide numerical 
examples to show that our algorithms are tractable for as many as three searchers 
and ten targets. For (GTP C ) we compare an algorithm based on fixed discretization 
schemes to an algorithm based on adaptive discretization schemes. Our results indi¬ 
cate that both methods produce equivalent solutions, but the adaptive discretization 
schemes require less computational time to acheive them. 

We compare three heuristic approaches that have the potential to be used 
as real-time solution methods onboard unmanned vehicles. We use two different 
polynomial based approaches, and one fixed-precision method to solve three randomly 
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generated problem instances. For these problem instances, the fixed-precision method 
offers the best overall performance. It consistently provides solutions with the best 
objective values and the lowest computational times. 

We also develop rate of convergence results for approximations to infinite di¬ 
mensional problems that involve a single discretization parameter for the space of 
decision variables and two discretization parameters for the objective function as a 
computational budget b tends to infinity. We note that the generalized optimal control 
problems we define in Chapter III are examples of these infinite dimensional problems. 
We show that superlinear and linear algorithms can achieve the same theoretical rate 
of convergence as that achieved by the “ideal” case of a finitely convergent algorithm. 
We also identify specific optimal discretization policies for both the superlinear and 
linear cases that achieve this best possible rate. Our analysis indicates that if a linear 
or superlinear optimization algorithm is used to solve the finite dimensional optimal 
control problem, with Euler’s method used to numerically solve the differential equa¬ 
tions and Simpson’s rule used to numerically approximate the spatial integration, 
then the best possible rate of convergence is b~ 2 ^. If a second-order Runge-Kutta 
method is used instead of Euler’s method to numerically solve the differential equa¬ 
tions, then the best possible rate of convergence is b~ l . If a fourth-order Runge-Kutta 
method is used instead of Euler’s method to numerically solve the differential equa¬ 
tions, then the best possible rate of convergence is 5 -4//3 . Finally, if an ideal method 
is used instead of Euler’s method to numerically solve the differential equations, then 
the best possible rate of convergence is b~ 2 . If “ideal” methods are used to solve the 
differential equations as well as evaluate the spatial integration, there is no benefit 
associated with increasing the discretization parameters N or M, so b = n, where n 
is the number of iterations of the optimization algorithm. The resulting asymptotic 
rates for a superlinear optimization algorithm with order 7 G (1, 00 ) and c G (0,1), 
and a linear optimization algorithm with rate constant c G (0,1) are c< b and c b , 
respectively. Based on our analysis, it appears that it is possible to improve the run- 
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time performance of our algorithms if a higher-order method such as Runge-Kutta or 
pseudo-spectra! is used in place of Euler’s method to numerically solve the differential 
equations. 

B. FUTURE WORK 

This dissertation suggests a number of possibilities for extension and future 
work. The first area we consider for additional work is the development of a real-time 
solution method. As discussed in Section V.B.3.d, the implementation of Algorithm 
V.6 used finite differences to estimate the required gradients. Future work should 
include explicit computation of all required derivatives, so that the implementation 
of Algorithm V.6 is on a more equal footing with those of the other two potential 
real-time methods. An opportunity for extension would be to use something other 
than algebraic polynomials as the basis for the searcher trajectories in Algorithms 
V.5 and V.6. One suggestion would be to use Chebyshev polynomials. In Caporale 
& Cerrato (2008) Chebyshev polynomials in conjunction with Chebyshev nodes were 
shown to provide excellent approximations to the solutions of linear partial differential 
equations, indicating that there might be potential for their use in Algorithms V.5 or 
V.6. Another suggestion would be to use Bezier curves in Algorithm V.6. Because 
Bezier curves are completely contained in the convex hull of their control points, they 
would be particularly well suited to an indirect curve fitting method where the control 
points could be determined based on the low level initialization as in Algorithm V.6. 

As discussed in Section VI.A above, there is potential for run-time improve¬ 
ment if higher order methods are used in place of Euler’s method to numerically solve 
the differential equations. The consistent approximation theory and implementable 
algorithms developed in Chapters III and V, respectively, could be extended based 
on using Runge-Kutta or pseudo-spectral methods instead of Euler’s method to nu¬ 
merically solve the differential equations, however, we foresee numerous technical 
challenges. 
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In this dissertation, we consider two different types of objective functions. An¬ 
other natural extension would be to consider additional types of objective functions. 
In the sections that follow, we present two alternative objective functions that, while 
potentially more difficult to solve, may be of operational interest in many situations. 


1. Minimize Expected Time Until First Detection 

We now consider a new performance metric for our optimization problem that 
is based on minimizing the expected time until the first detection of all the targets. 
In Chapter II, we derived (II.3), which is an expression for the probability that the 
k th searcher does not detect the I th target during [0, t], t e [0,1] We define the random 
variable T as the time of first detection of the I th target, and note that 

P ({ T > t}) = q k ’ l (t, a) = exp (- jT r k - 1 {x k (s), y l (s; a)) ds^j = 1 - F(t), (VI.l) 


where F(t) is the cumulative distribution function for the random variable T. Then 
for the case of one searcher and one target the conditional expected time of first 
detection given a particular target trajectory is given by 


E[T\a\ = / tf T (t)dt, 


(VI.2) 


where frit) is the probability density function of the random variable T. We note 
that the support of the random variable T is [0,1], so E[T\a\ is finite. We also note 
that we can use (VI.l) and integration by parts to show that (VI.2) can be re-written 
in terms of q k,l (t ; a) as follows 

E[T\a] = j q k,l (t-,a)dt = J exp J' r k ’ 1 (x k (s),y l (s;a)) ds'j dt 

= [ p({T>t})dt = P({T>t})t\ 1 - f t ( -hit )) dt 




(VI-3) 


If we assume that the searchers make independent detection attempts, then for the 
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case of multiple searchers the conditional expected time of first detection given a 
particular target trajectory is given by 



If we assume that the random variables that the target motion is conditioned upon 
are independent across targets, then the expected time of first detection of the I th 
target is given by 


ds J dt 

2. Herding Formulation 

Throughout this dissertation, we have focused on detecting the potential threats 
to the HVU. Once the potential threats are located, the searchers could become de¬ 
fenders who now select controls such that the potential attackers are herded away 
from the HVU. In a herding model, we assume a different target motion model where 
the targets’ motion is defined by coupled dynamics, which we describe in detail below. 
Then, a natural metric for the degree of success in herding would be to maximize the 
minimum distance between the HVU and any attacker at any instance of time in the 
planning horizon. This herding success metric is given by 


(j) 1 (a 1 ) da 1 . (VI.7) 



(x°(t) -y[ (t;a l )) 2 + (x° 2 {t) - y\ (t;a l )) Z 4> l [a 1 ) da 1 \ . (VI.8) 


m 2 l 


max < nun / 

z fe O l M Ja l &A L 

In order to formulate the dynamics between the defenders and the attackers, 
we first define some necessary quantities. First, we define 


ran 9 ^ d {t) = yj(x’ftt) - y\ (t; a 1 )) 2 + (x k 2 [t) - y l 2 (t; a 1 )) 2 , (VI.9) 
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which gives the range between the I th attacker and the k th defender at time t, given 
a particular attacker trajectory. We also define 


range l HVU {t) = \J (x?(t) - y[ (t; a 1 )) 2 + (x$(t) - y l 2 (t; a 1 )) 2 , (VI. 10) 


which gives the range between the I th attacker and the HVU at time t, given a 
particular attacker trajectory. Next, we define 


9 ad^) — tan 1 


Aif) -y l 2 
(t) -y[ (t; a 1 ) J 


(VI-11) 


which gives the angle between the I th attacker and the k th defender at time t, given 
a particular attacker trajectory. Finally, we define 


9 Hvu(t) — t an 


S°(*) -y l 2 

xi(t) ~y\ i^ a ) J 


(VI.12) 


which gives the angle between the I th attacker and the HVU at time t, given a par¬ 
ticular attacker trajectory. Then, the dynamics of the I th attacker given a particular 
realization of the random variable a are given by 

I< 

y[(t) = 


COS d l HV uit) 


k= 1 
K 


y l 2 (t) = 


™nge l a k d (t) 

sin ^W 


+ wl 


[range l HVU {t)\ 

sin 9 l HV u(t) 


k= i \range l £(t)] [™nge l HVU {t)\ 


2 ’ 


(VI.13) 


(VI. 14) 


where w\ and w l 2 are user specified weights that define the relative importance of the 
two components of the attacker dynamics. The component associated with w[ is due 
to the herding effects of the defenders. In a manner similar to the approach used in 
Lu (2006), we assume that the I th attacker moves in a straight line away (hence the 
minus sign in front of w \) from the weighted sum of “influences” of the k defenders, 
and that the I th attacker’s velocity is bounded by the inverse of the distance between 
it and the k th defender squared. The component associated with w l 2 is due to the 
attacker’s desire to get close to the HVU in order to facilitate his attack. Again, we 
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assume that the I th attacker moves in a straight line towards (hence the plus sign in 
front of w l 2 ) the HVU, and that the I th attacker’s velocity is bounded by the inverse 
of the distance between it and the HVU squared. 

One way that w[ and w l 2 could be specified would be to define a distance called 
panic 1 that gives the range from the HVU that the I th attacker begins to worry about 
being detected by the defenders. In this case w\ and w l 2 would be given by 

, panic 1 

w i = - 1 - 

range HVU 

„j _ range l HVU 

^2 • / 
panic 1. 


(VI.15) 
(VI.16) 
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VII. 


APPENDIX: MATHEMATICAL 
BACKGROUND 


Throughout this dissertation, we consider the Lipschitz continuity of functions 
relative to the set H, and the differentiability of functions on the set H°, relative to 
the set H as defined in Polak (1997). For the sake of completeness, we include the 
pertinent definitions from pp. 652-656 of Polak (1997). We begin with the definitions 
related to continuity. 

Definition VII. 1. Let V be a real nornred space and let S be a convex subset of V. 

(i) A function / : V —> M m is said to be continuous at a point x G V, if, for every 
5 > 0, there exists a p > 0 such that 

\\f{x')^f{x)\\<8yx'eB{x,p) r (VII.1) 

o 

where B(x, p) = {x' G M n j \\x' — x|| < p}. A function / : V —> is said to be 

continuous (continuous on S ) if it is continuous at all x G V (x E S). 

(ii) A function / : V —> M m is said to be continuous, relative to S (S'-continuous), 
if, for every x G S and for every <5 > 0, there exists a p > 0 such that 

11/00 - f{x) II < S,W G B(x,p) n S. (VII. 2) 


Next we state definitions related to differentiability. 


□ 


Definition VII.2. Let V be a real normed space. 


(i) We will say that a continuous function / : V —> M m is Gateaux differentiable at 
a point x* G V, if there exists a bounded linear operator f x {x*) '■ V —> M m such 
that, for every 5x G V, 


lim || fjx* + Ahx) - f(x*) - A f x {x*)Sx 
A4-0 A 


(VII.3) 


We will call f x (x*) the Gateaux derivative of /(•) at x*. 

We will say that a continuous function / : V —> M m is Gateaux differentiable on 
a subset S of V, if it is Gateaux differentiable at all x G S. 
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(ii) We will say that a continuous function / : V —> M m is Frechet differentiable at 
a point x* G V, if it is Gateaux differentiable at x*, with Gateaux derivative 
f x (x*), and, in addition, the Gateaux derivative has the property that 


lim 

Sx—>0 


|| f(x* + Sx) - f(x*) 

IIM 


fx(x*){Sx) II 


0 . 


(VII.4) 


In this case, we will also call f x (x*) the Frechet derivative of /(•) at x*. 

We will say that a continuous function / : V —>• M m is Frechet differentiable on 
a subset S of V, if it is Frechet differentiable at all x € S. 


(iii) Let S' C 5 be two convex subsets of V. We will say that a continuous function 
/ : S —> M m is Gateaux differentiable, relative to S, (Gateaux S-differentiable) 
at a point x* G S', if there exists a bounded linear operator f x (x*) : V —>• M m 
such that, for every Sx G V such that x* + A *Sx G S, for some A* > 0, 


lim II f(x* + XSx) - f(x*) - A f x (x*)Sx 
A4-0 A 


(VII.5) 


We will call f x (x*) the Gateaux ,S'-derivative of /(•) at x*. 

We will say that /(•) is Gateaux S'-differentiable on S', if it is Gateaux S- 
differentiable at all igS'. 


(iv) Let S' C S be two convex subsets of V. We will say that a continuous function 
/ : S —> is Frechet differentiable, relative to S, (S'-differentiable) at a point 
x* G S', if it is Gateaux S’-differentiable at x*, with Gateaux S-derivative f x (x*), 
and, in addition, the Gateaux S-derivative has the property that 


lim 

Sx —>-0 
x*+Sx€S 


|| f(x* + 5x) - f(x*) 

IHI 


f x (x*){Sx) || 


0 . 


(VII.6) 


In this case, we will also call f x (x*) the Frechet S-derivative of /(•) at x*. 


We will say that /(•) is S-differentiable on S', if it is Frechet S-differentiable at 
all a: G S'. 


□ 

The next two paragraphs are a restatement of the beginning of Section 5.6.2 
from Polak (1997), but are included here for the sake of completeness to provide a 
detailed explanation of what it means for a function to be continuous relative to the 
set H and Gateaux differentiable at rj G H 11 , relative to the set H. 
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Based on Definition VII.1, we can now give a detailed explanation of what it 
means for a function to be continuous relative to the set H. Let S be a convex subset 
of H. Then, according to Definition VII. 1, a function / : S —» M n is continuous 
relative to H. at a point rj G S, if, for any sequence with r/, G H for all 

i G N, such that rj, —>• rj, as i —> oo, /(r^) —* /(/?), as i —* oo. It should be noted, 
however, that if {? 7 j}? 2 0 is any arbitrary sequence converging to rj, then it is possible 
that r/j H for all i G N. Then, we cannot claim that f{j]j) —> /(r?), as i —>■ oo, 
and hence we cannot necessarily conclude that /(•) is continuous at rj. To keep our 
terminology as concise as possible, as in Definition VII. 1, given S, a convex subset of 
H. we will say that a function / : S —» M n is H-continuous when we mean that it is 
continuous relative to H. 

Similarly, based on Definition VII.2 we can give a detailed explanation of what 
it means for a function to be Gateaux differentiable at r) G H°, relative to the set H. 
We begin by noting that given any r] G H° and any dr] G H 00,2, there always exists 
a A > 0 such that rj + A 5r) G H. Then it follows from Definition VII.2(iii) that a 
function / : H — > M n is Gateaux differentiable at rj G H°, relative to H. if there exists 
a continuous linear map Df(r /; •) from H 00j 2 into M n , such that for all Srj G LGo,2, 


Um f(v + A^?y) - f(rj) - A Dfjrf, 5rj) 
A4-0 A 


(VII.7) 


As a result, if / : H —y M n is Gateaux differentiable at r/ G H°, relative to H. it is also 
Gateaux differentiable at rj. It also follows from Definition VII.2(iv) that / : H —>• 
is Frechet differentiable relative to H (Frechet H-differentiable) at rj G H n , if it is 
Gateaux differentiable at rj and the Gateaux differential Df(r]; •) has the property 


II fW) - f(v) ~ Df (rj; rj - 77)|| 

,.!-p -VWG- = 0 ' 


(VII.8) 


Then by definition, a function / : H —y M n is Gateaux/Frechet H-differentiable on 
H°, if it is Gateaux/Frechet H-differentiable at every rj G H n . 
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